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

HAL
open science

# The influence of gas pore pressure in dense granular flows: numerical simulations versus experiments and implications for pyroclastic density currents

Alvaro Aravena, Laurent Chupin, Thierry Dubois, Olivier Roche

► To cite this version:

Alvaro Aravena, Laurent Chupin, Thierry Dubois, Olivier Roche. The influence of gas pore pressure in dense granular flows: numerical simulations versus experiments and implications for pyroclastic density currents. Bulletin of Volcanology, 2021, 83 (11), 10.1007/s00445-021-01507-7. http://dx.doi.org/10.1007/s00445-021-01507-7

HAL Id: hal-03408294

https://uca.hal.science/hal-03408294

Submitted on 29 Oct 2021

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1

# The influence of gas pore pressure in dense granular flows: numerical simulations versus experiments and implications for pyroclastic density currents

2

3

Alvaro Aravena¹,², Laurent Chupin², Thierry Dubois², Olivier Roche¹

4

¹Laboratoire Magmas et Volcans, Université Clermont Auvergne, CNRS, IRD, OPGC, Clermont-Ferrand, France.

5

6

²Laboratoire de Mathématiques Blaise Pascal, Université Clermont Auvergne, CNRS, Clermont-Ferrand, France.

7

8

Corresponding author: Alvaro Aravena (alvaro.aravena@uca.fr).

9

## Abstract

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

We investigate the influence of gas pore pressure in granular flows through numerical simulations on horizontal and low-angle inclined surfaces. We present a two-phase formulation that allows description of dam-break experiments considering high-aspect-ratio collapsing columns and depth-dependent variations of flow properties. The model is confirmed by comparing its results with data of analogue experiments. The results suggest that a constant, effective pore pressure diffusion coefficient can be determined in order to reproduce reasonably well the dynamics of the studied dam-break experiments, with values of the diffusion coefficient consistent with experimental estimates from defluidizing static columns. The discrepancies between simulations performed using different effective pore pressure diffusion coefficients are mainly observed during the early acceleration stage, while the final deceleration rate, once pore pressure has been dissipated, is similar in all the studied numerical experiments.

However, these short-lasting discrepancies in the acceleration stage can be manifested in large differences in the resulting run-out distance. We also analyze the pore pressure at different distances along the channel. Although our model is not able to simulate the under-pressure phase generated by the sliding head of the flows in experiments and measured beneath the flow-substrate interface, the spatio-temporal characteristics of the subsequent over-pressure phase
are compatible with experimental data. Additionally, we studied the deposition dynamics of the granular material, showing that the timescale of deposition is much smaller than that of the granular flow, while the time of the deposition onset varies as a function of the distance from the reservoir, being strongly controlled by the surface slope angle. The simulations reveal that an increment of the surface slope angle from 0° to 10° is able to increase significantly the flow run-out distance (by a factor between 2.05 and 2.25, depending on the fluidization conditions). This has major implications for pyroclastic density currents, which typically propagate at such gentle slope angles.

34

# 1. Introduction

Pyroclastic density currents (PDCs) are gravity-driven flows of hot particles (pyroclasts and lithic fragments) and gas (Druitt 1998; Branney and Kokelaar 2002; Dufek et al. 2015; Dufek 2016) generated by the partial or total collapse of an eruptive column or a volcanic dome. They exhibit a wide range of particle concentration, temperature and grain size distribution, and two physical regimes can be recognized as end-members (Branney and Kokelaar 2002; Burgisser and Bergantz 2002; Dufek 2016): dilute and dense flows, which may occur alone (e.g. dilute turbulent flow) or simultaneously (e.g. dense base and overriding turbulent part) during the propagation of a PDC. The dilute component of PDCs consists in a turbulent suspension with a solid concentration of the order of 1 vol. % or less dominated by the interaction between solid particles and the interstitial gas, while the dynamics of the dense component of PDCs, if present, is typically dominated by particle-particle interaction and by friction with the topography, presenting a solid concentration of the order of 30 – 60 vol. % (Lube et al. 2020). A transport regime of PDCs characterized by clusters at intermediate particle concentrations (i.e., a few vol.% to ~30 vol.%) has been recognized recently and may be present for instance in a transitional zone between a dense base and an upper dilute turbulent part (Breard et al. 2016;
Fullmer and Hrenya 2017; Lube et al. 2020). Because of their high propagation velocities, dynamic pressures and temperatures, PDCs can devastate urbanized zones, being one of the most hazardous processes associated with volcanic eruptions (Druitt 1998; Branney and Kokelaar 2002; Cole et al. 2015; Neri et al. 2015). Thus, deciphering the factors controlling their dynamics and the expected run-out distance is of paramount importance for volcanic hazard assessment. Although much attention has been paid to the study of the long run-out distance that characterize some PDCs (Bursik and Woods 1996; Branney and Kokelaar 2002; Kelfoun 2011; Roche et al. 2016, 2021; Shimizu et al. 2019; Giordano and Cas 2021), several aspects remain poorly understood. PDC run-out distance is the result of a series of concomitant processes whose relative efficiency is influenced by the flow properties (e.g. solid particle concentration, volume, speed and temperature) and the regional slope (Valentine et al. 2011), and include: the interaction with the surrounding atmosphere (e.g., air entrainment and heat transfer; Benage et al. 2016), the rheological effect of interstitial pore fluid pressure (Druitt et al. 2007; Roche 2012), and the interplay between the flow base and the substrate, where different processes may occur, such as erosion (Cas et al. 2011; Bernard et al. 2014; Farin et al. 2014), self-channelization (Brand et al. 2014; Gase et al. 2017), self-fluidization (Breard et al. 2018; Chédeville and Roche 2018), and pyroclast deposition (Branney and Kokelaar 2002).

In particular, within the pyroclastic mixture, and especially at the impact zone of a collapsing fountain (Sweeney and Valentine 2017; Valentine and Sweeney 2018; Valentine 2020; Fries et al. 2021), the differential motion between the interstitial gas (flowing relatively upwards) and the solid particles (moving relatively downward) is able to generate pore pressure, which counterbalances the weight of the particles, reduces friction and thus increases run-out distance (Iverson 1997; Savage and Iverson 2003; Goren et al. 2010; Roche 2012; Rowley et al. 2014; Breard et al. 2019a). The temporal evolution of pore pressure, and thus its effective influence on run-out distance, depends on the balance between some source mechanisms (e.g. gas
ingestion, differential gas-particle motion caused by particle settling) and pore pressure
diffusion, which is in turn controlled by the properties of the PDC material. In fact, slow gas
pressure diffusion is favored by thick pyroclastic flows and by grain size distributions
dominated by fine particles that confer low hydraulic permeability (Druitt et al. 2007; Burgisser
2012; Roche 2012; Breard et al. 2019b).

In this work we address the influence of pore pressure on the propagation of granular flows
through numerical simulations. In particular, we present a two-phase model, built on the
formulation presented by Chupin et al. (2021), which accounts for the effect of pore pressure
on the dynamics of granular flows and allows us to simulate collapsing columns in the dam-
break configuration and the subsequent flow propagation on horizontal and low-angle inclined
surfaces. The column height and aspect ratio adopted in our numerical simulations (40 cm and
2, respectively) were selected to allow model confirmation by comparing numerical results with
published experimental data (cf. Valentine 2019; Esposti Ongaro et al. 2020) of collapsing
columns over a horizontal surface (Roche et al. 2010). Note that we use the term *confirmation*
instead of *validation* following the framework presented by Esposti Ongaro et al. (2020).
Numerical results also allow us to explore some key physical aspects controlling the dynamics
of granular flows (e.g. pore pressure spatio-temporal evolution and flow deposition), which are
often difficult to measure in time across the entire spatial domain of analogue experiments.
Moreover, adopting a set of input conditions calibrated using experimental data, we performed
additional simulations considering collapsing columns on low-angle inclined rigid surfaces, in
order to test the coupled effect of pore pressure and topography on the propagation of granular
flows. Compared to previous efforts to address numerically the influence of pore pressure in
the propagation dynamics of PDCs (Gueugneau et al. 2017), which are based on depth averaged
models, our model has some relevant strengths: it allows us to study high-aspect ratio collapsing
columns and to describe depth-dependent variations of the flow properties.
This article consists of five sections. In Section 2 we describe the experimental configuration considered in this paper. In Section 3, we present the numerical model adopted (Section 3.1), its confirmation by comparing numerical results with those of analogue experiments (Section 3.2), and then we describe the results of simulations performed considering low-angle inclined surfaces (Section 3.3). Finally, in Sections 4 and 5 we present the discussion and concluding remarks of this article.

107

## 2. Experimental configuration

In order to test and confirm our model, we considered the experimental data presented by Roche et al. (2010). The benchmark experiment consists in the release of a fluidized granular column into a horizontal, smooth channel (note that the term *fluidization* is used here to refer to the presence of a vertical flow of air able to counterbalance the bed weight, and it is not related to the presence of other fluid phases such as water). The dynamics of the dam-break experiment, which was measured using high-speed cameras and pressure sensors located at different positions along the horizontal channel, can be decomposed into three stages: (1) a quick phase of initial acceleration, (2) propagation of the front at nearly constant velocity, and (3) deceleration of the flow and front stopping. Roche et al. (2010) and Chupin et al. (2021) also pointed out a final stage of very slow propagation of granular material in the flow body after the front stopped. The experimental apparatus includes a reservoir of 20 cm length and 10 cm width, and a channel of 3 m length and 10 cm width. Initially, the particles are introduced into the reservoir (column height of 40 cm) where an air flow is supplied from below in order to generate fluidization and the related pore pressure. This simple configuration aims to mimic particle-gas differential motion generated through various means, including particle settling (Chédeville and Roche 2018; Breard et al. 2018; Valentine and Sweeney 2018; Fries et al. 2021). Roche et al. (2010) tested two fluidization conditions by adjusting the supplied air
velocity: (1) imposing the minimum fluidization velocity ($U_{mf}$) or (2) imposing the minimum bubbling velocity ($U_{mb}$). The minimum fluidization velocity $U_{mf}$ (~0.8 cm/s in the experiments; Roche et al. 2010) guarantees that the bed weight is counterbalanced by the drag of the interstitial air flow on the particles, and the granular bed is not expanded. On the other hand, at $U_{mb}$ (~1.3 cm/s in the experiments; Roche et al. 2010), the bed weight is counterbalanced and the granular network is expanded. In order to trigger column collapse, at $t = 0$, a sluice gate is opened rapidly (< ~0.1 s), allowing to release the granular material, which propagates laterally along the horizontal channel during about 1.3 s. As our numerical model treats incompressible flows, we compare our results with the analogue experiment performed using the minimum fluidization velocity ($U_{mf}$), that is, when the bed is not expanded. The particles used in these experiments were spherical glass beads with a grain size range of 60-90 µm (monodisperse size distribution, mean of 75 µm) and a density of $\rho_s$ = 2500 kg/m³. Note that more complex particle shapes and size distributions are able to control pore pressure diffusion in granular flows by affecting porosity and mixture permeability (Wilson 1984; Burgisser 2012; Breard et al. 2019b and references therein), and that Breard et al. (2019b) showed that the effective particle size regarding fluidization and pore pressure diffusion is the Sauter diameter, which is very close to the mean diameter for subspherical particles such as we considered. The resulting granular column had a bulk density of $\rho_b$ = 1450 ± 50 kg/m³ (i.e. pore volume fraction of $\epsilon$ = 0.42 ± 0.02). Additionally, we can calculate the theoretical hydraulic diffusion coefficient $\kappa_T = k/(\epsilon\mu\beta)$, where k is hydraulic permeability, $\mu$ is gas dynamic viscosity and $\beta$ is gas compressibility. In case of a perfect gas, $\beta = 1/P_i$, where $P_i$ is the initial pore pressure, which is about equal to the atmospheric pressure. Considering that $k \sim 10^{-11}$ m², $\epsilon \sim 0.42$, and $\mu \sim 1.8 \times 10^{-5}$ Pa s, we obtain $\kappa_T \sim 0.13$ m²/s. However, it is worth noting that this value is one order of magnitude larger than the estimates of diffusion coefficient given by Roche et al. (2010) ($\kappa \sim 0.01$ m²/s), which are based on
experimental measurements on static defluidizing beds and are shown to increase with the bed height. The reason explaining this discrepancy is unknown and is discussed below.

## 3. Numerical simulations

### 3.1 Mathematical modelling and numerical schemes

Based on the numerical model presented by Chupin et al. (2021), we constructed a new model able to consider the effect of pore pressure and reproduce the experimental configuration adopted by Roche et al. (2010). We consider the collapse of a granular mass over a planar rigid surface with inclination angle $θ$ varying from horizontal up to 10°. As the laboratory experiments have been performed in a narrow channel (10 cm wide and 3 m long; Roche et al. 2010), we consider the problem as mainly two-dimensional. Note that we neglect the effects of the lateral walls.

The granular medium, which is a mixture of air and glass beads, is described by an incompressible flow with a $μ(I)$-rheology (Jop et al. 2006). In this rheological model, which has been widely adopted to describe dense granular flows (Gray and Edwards 2014; Ionescu et al. 2015), the dynamics of the granular flow is governed by the mass and momentum conservation laws

$$ \rho(\partial_t u + u \cdot \nabla u) = \text{div} T + \rho g, \qquad (1) $$

$$ \text{div} u = 0, \qquad (2) $$

$$ \partial_t \rho + u \cdot \nabla \rho = 0, \qquad (3) $$

where $u$ is the material velocity, $g$ is an external force (gravity), $T$ is the total stress tensor, and $ρ = φρ_s$ is the bulk density, where $ρ_s$ and $φ$ are the particle density and average volume fraction, respectively. In our simulations, based on the experimental data described by Roche et al. (2010), we use $ρ_s = 2500 \text{ kg/m}^3$ and $φ = 1 - ε = 0.58$. The pressure is given by $p =$
173 $-\frac{1}{3}\text{tr } \mathbf{T}$, so that the deviatoric stress $\mathbf{T}'$ ($\mathbf{T} = -p\mathbf{Id} + \mathbf{T}'$) has to be prescribed in order to close equations (1) – (3).

175 Modeling the granular flow with equations (1) – (3) entails that the total pressure $p$ is the sum of the solid (effective) pressure $p_s$, due to force chains of glass beads, and the pore pressure $p_f$, 176 due to the presence of air between particles. Therefore, in order to account for the presence of 177 air between glass beads, the pore pressure and its effect on granular flows should be modeled. 178

179 Following Iverson and Denlinger (2001), the pore pressure diffuses and is advected with the granular mass so that $p_f$ satisfies the balance equation

$$ \partial_t p_f - \text{div}(\kappa \nabla p_f) + \mathbf{u} \cdot \nabla p_f = 0, \quad (4) $$

181 where $\kappa$ is the diffusion coefficient. The knowledge of $p_f$ through equation (4) permits us to 182 define the effective pressure as $p_s = p - p_f$.

183 In the $\mu(I)$-rheology (Jop et al. 2006), the deviatoric stress $\mathbf{T}'$ is given by

$$ \mathbf{T}' = \mu(I)p_s \frac{\mathbf{D}(\mathbf{u})}{|\mathbf{D}(\mathbf{u})|}, \quad (5) $$

184 where $\mathbf{D}(\mathbf{u}) = \frac{1}{2}(\nabla\mathbf{u} + \nabla\mathbf{u}^t)$ is the strain rate tensor and $|\mathbf{D}(\mathbf{u})|^2 = \frac{1}{2}\sum_{i,j} \mathbf{D}(\mathbf{u})_{i,j}^2$. The friction coefficient $\mu(I)$ depends on the inertial number $I$, namely

$$ \mu(I) = \mu_s + \frac{\mu_\infty - \mu_s}{1 + I_0/I} \quad \text{and} \quad I = \frac{2d|\mathbf{D}(\mathbf{u})|}{\sqrt{p_s/\rho_s}}. \quad (6) $$

187 In equation (6), $d$ is the particle diameter, $I_0$ is a dimensionless number, $\mu_s = \tan(\alpha)$ with $\alpha$ representing the static internal friction angle of the granular material, and $\mu_\infty \ge \tan(\alpha)$ is an 188 asymptotic value of the friction coefficient for large inertial numbers. By combining equation 189 (6) and equation (5) (see Chupin et al. (2021) for the details), we can rewrite the expression of 190 the tensor $\mathbf{T}'$ in regions where $\mathbf{D}(\mathbf{u}) \ne 0$ as

$$ \mathbf{T}' = 2\eta(|\mathbf{D}(\mathbf{u})|, p_s)\mathbf{D}(\mathbf{u}) + \tan(\alpha) p_s \frac{\mathbf{D}(\mathbf{u})}{|\mathbf{D}(\mathbf{u})|}, \quad (7) $$
$$ \eta(|\mathbf{D}(\mathbf{u})|, p_s) = \frac{(\mu_\infty - \tan(\alpha))p_s}{2|\mathbf{D}(\mathbf{u})| + \frac{l_0}{d}\sqrt{p_s/\rho_s}}. \quad (8) $$

With this formulation, the μ(I)-rheology appears to be a viscoplastic rheology with a Drucker-Prager plasticity criterion (see also Jop et al. 2006 and Ionescu et al. 2015) and a spatio-temporal variable viscosity, which is a fundamental aspect of our study. In order to treat the non-differentiable definition of the tensor T', due to the presence of the 1/|D(u)| term that is singular in the absence of strain rate, a projection scheme is applied (Chalayer et al. 2018; Chupin et al. 2021). The projection procedure avoids a need to resort to any regularization technique and allows to accurately capture the rigid zones, i.e., the regions where no deformation occurs.

As in Chupin et al. (2021), the presence of the ambient gas (i.e. the air outside the flow) is taken into account. The granular flow and the ambient air flow are separated by an interface transported by the velocity field. A level-set function Φ (see Osher and Fedkiw (2001) for instance), initially defined as the signed distance to the interface, is used to describe the limit between the granular flow and the ambient gas. The computational domain is split so that Φ < 0 corresponds to the granular flow, Φ > 0 the ambient gas and Φ = 0 the interface. The level-set function satisfies the equation

$$ \partial_t \Phi + \mathbf{u} \cdot \nabla \Phi = 0. \quad (9) $$

The ambient flow (Φ > 0) is also governed by equations (1) – (3) but with a Newtonian rheology, namely T' = 2η_f D(u) where η_f is the air dynamic viscosity, and a mass density ρ = ρ_f. Note that the pore pressure equation (4) has a meaning only inside the granular flow, that is where Φ < 0. In order to solve an equation valid over the whole computational domain, the diffusion coefficient κ takes a very large value (≈ 10¹⁶ m²/s) outside the granular flow so that p_f is extended to zero outside of the granular flow.
Coulomb friction boundary conditions are applied on the vertical backwall of the reservoir and on the bottom of the channel on which the granular medium slides, that is

$$u_n = 0 \text{ and } T_t = -\eta_b u_t - \tan(\alpha_b) (-T_n)^+ \frac{u_t}{|u_t|}, \quad (10)$$

where $u_n = u \cdot n$ ($n$ being the unit outward normal vector to the domain boundary) is the normal velocity and $u_t = u - u_n n$ the tangential one. Similarly, for the stress we have $T_n = (T \cdot n) \cdot n$ and $T_t = T \cdot n - T_n n$. In all simulations reported in this paper, the friction angle on the vertical backwall of the reservoir and on the bottom of the channel ($\alpha_b$) was set to 15° (Chupin et al. 2021).

At time $t = 0$, the pore pressure $p_f$ is initialized in the reservoir as 90% of the weight of the particles, that is, $p_f(x, y) = 0.9\rho g(H - y)$ ($H$ being the height of the initial column), which agrees with experimental measurements (Montserrat et al. 2012). Equation (4) is supplemented with Neumann boundary conditions. On the bottom of the channel inside the reservoir, that is for $x \in [-20 \text{ cm}, 0 \text{ cm}]$ and $y = 0$, the constant air flux imposed in the experiment is modeled with a constant pressure gradient $\frac{\partial p_f}{\partial n} = -0.9\rho g$. Everywhere else on the domain boundary, a homogeneous Neumann boundary condition $\frac{\partial p_f}{\partial n} = 0$ is applied.

Equations (1)–(3) and (4) are discretized in space with second order finite volume schemes on a staggered grid. A bi-projection scheme (Chalayer et al. 2018) is applied for the temporal discretization. The level-set transport equation (9) is solved with a RK3 (third-order Runge-Kutta) TVD (total variation diminishing) scheme coupled with a fifth-order WENO (weighted essentially non-oscillatory) scheme. We also apply a reinitializing algorithm in order to maintain the level-set function to the signed-distance function to the interface between the granular and the ambient flows. Details on the numerical schemes are provided in Chupin and Dubois (2016), Chalayer et al. (2018), and Chupin et al. (2021). In order to remedy to a lack of resolution of the small scale structures outside the granular flow, due to the small value of the
237

air viscosity, a subgrid scale (Smagorinsky) model is used (Smagorinsky 1963). This results in
enhancing the viscosity of air $\eta_f$ by adding a local eddy viscosity defined as $(C_s h)^2 |D(\mathbf{u})|$, where $h$ is the mesh size and $C_s \in [0.1, 0.2]$ is the Smagorinsky constant. We used hereafter the value 0.1 for $C_s$ in all simulations.

238

239

240

241

The code, written in F90, is parallel: the PETSc library (Balay et al. 2018a, 2018b) is used to solve the linear systems and to manage data on structured grids while communications between processes are explicitly written with MPI subroutines.

242

243

244

**3.2 Simulations on horizontal planes**

245

246

247

248

249

250

251

252

253

As a first step, we performed a set of simulations on horizontal surfaces using the model described in Section 3.1 and considering different values of the pore pressure diffusion coefficient ($\kappa$). Here we focus on the results of simulations that agree reasonably well with the experimental results presented by Roche et al. (2010) in terms of run-out distance, temporal evolution of the front position and profile of the flow surface, that is, with pore pressure diffusion coefficients ranging from 0.015 to 0.035 m²/s. Interestingly, these values agree with those determined experimentally by measuring the timescale of pore pressure diffusion in static columns of heights of ~15-35 cm. For comparison purposes, we also include the results of a test simulation performed using the theoretical hydraulic diffusion coefficient ($\kappa_T$~0.13 m²/s).

254

255

256

257

258

259

Note that we consider constant effective diffusion coefficients (with values of 0.015, 0.025, and 0.035 m²/s), while the diffusion coefficient likely varies during the analogue experiment due to granular material dilation and compaction. To reproduce the initial setup of the benchmark analogue experiment, in our simulations the initial height of the collapsing column is 40 cm and the initial width is 20 cm. The height of the computational domain is 45 cm, adopting a grid with 128 cells in the vertical direction (i.e. cell size of 3.5 mm).

260

261

Our model tends to under-estimate the deposit thickness in proximal domains (from <5% up to ~25%) and to over-estimate the deposit thickness in distal domains (Figs. 1a-c and 2). Still, the
general shapes of the simulated final profiles of the deposits are very similar to that of the benchmark analogue experiment, i.e. profiles dipping gently downstream and with the maximum thickness located in the channel near the reservoir. This deposit shape differs clearly from that of non-fluidized granular flows, whose maximum thickness is located at the back-wall of the reservoir while the thickness decreases monotonically with distance (Roche et al. 2010; Ionescu et al. 2015). Moreover, numerical results reproduce reasonably well the three phases of propagation described by Roche et al. (2010) (Fig. 3), and the relative duration of each phase as well as the flow duration are consistent with the benchmark experiment. The dynamics of gate opening in the analogue experiment slightly affects flow propagation during the initial acceleration phase, which may explain the differences in the initial front velocity during about 10% of the simulation duration (Figs. 3). An effective diffusion coefficient ($\kappa$) of about 0.015 m$^2$/s reproduces the experimental run-out distance, whereas a larger values of diffusion coefficient reproduce better the maximum thickness of the deposit. Note that, because our simulations are not able to describe flow thicknesses lower than 3.5 mm (i.e. the cell size used in numerical simulations), we compare our results with filtered experimental data, that is, with no consideration of flow thicknesses below this threshold (see Fig. 3 and its caption). On the other hand, the use of the theoretical value of the diffusion coefficient (i.e. 0.130 m$^2$/s) fails completely in reproducing the propagation dynamics of the benchmark experiment, under-estimating significantly the run-out distance (<60% of the run-out distance measured in the benchmark experiment; Figs. 1d, 2 and 3).

The simulations performed using $\kappa$ = 0.015 and $\kappa$ = 0.035 m$^2$/s give rise to differences of about 15% in the maximum velocity reached by the flow front (Fig. 4). The phase of velocity increase lasts ~17-22% of the whole propagation time, while the constant-velocity stage, which is slightly longer for simulations with low diffusion coefficients, represents ~15-25% of the total propagation time. Most of the propagation time of the granular flows (about 60-70%) is
associated with the phases of deceleration and front stopping. Our results of maximum velocity ($u/\sqrt{gH} \sim 1.0 - 1.15$) are consistent with the results presented by Roche et al. (2010), which further confirm the validity of our model. Note that the initial phase is characterized by the same acceleration in all the simulations, and the main differences between our simulations are observed in the absolute duration of this stage (and thus in the velocity reached by the flow front, Fig. 4). Another interesting result is that the velocity decrease during the final phase occurs at a similar rate in all the simulations (deceleration of ~0.23g). This shows that the differences in the front velocity during early phases of flow propagation are the cause of the different run-out distances, while negligible differences are observed in the dynamics of the deceleration stage. This is consistent with the fact that, once the initial pore pressure is completely dissipated, granular flows have the same rheological behaviour. On the other hand, as observed by Roche (2012), during the stopping stage the run-out distance increases with time to the power of 1/3.

The evolution of the pore pressure (Fig. 5) is the result of the coupled effect of diffusion, which occurs at a rate controlled by $\kappa$, and advection, controlled by flow velocity and thus in turn influenced by $\kappa$. Basal pore pressure undergoes an initial phase dominated by advection (i.e. advance of the iso-pressure fronts, which indicate the position along the x-axis at which specific values of pore pressure are reached as a function of time; Fig. 6), and a later phase dominated by diffusion (i.e. recession of the iso-pressure fronts) until stationary conditions are reached (Fig. 6). In particular, the simulations with $\kappa = 0.025 \text{ m}^2/\text{s}$ and $\kappa = 0.035 \text{ m}^2/\text{s}$ show a smooth, gradual transition between both phases. Instead, some of the iso-pressure fronts (Fig. 6) for the simulation with $\kappa = 0.015 \text{ m}^2/\text{s}$ show a significantly longer phase dominated by advection and then an abrupt decrease of pore pressure near the front. This rapid pore pressure decrease is favored by the small thickness of the granular flow at the front, while in the other cases such small flow thicknesses are reached while the flow is already defluidized.
In order to further compare our results with experimental data (Roche et al. 2010), we studied the pore pressure signal at different points along the channel base. Note that the under-pressure phase measured in experiments beneath the flow-substrate interface during the passage of the sliding head of the flow cannot be computed in our numerical simulations. Still, in Figure 7a-c we show the evolution of the modeled basal pore pressure at specific points along the x-axis, and we also display the differential pressure measured in the benchmark experiment (in the channel base at x = 0.2 m). The experimental data show that the passage of the flow front at a given point is followed by a short under-pressure phase and a later and longer over-pressure phase. Roche et al. (2010) propose that the under-pressure stage is mainly caused by the basal slip boundary condition and possibly by dilatancy processes (Garres-Díaz et al. 2020; Bouchut et al. 2021), which is supported by simulations (Breard et al. 2019a). Moreover, the minimum value reached during the under-pressure phase was empirically correlated to the slip velocity ($u_{slip}$; Roche 2012). The over-pressure phase would be instead dominated by compaction and advection of pore pressure within the granular flow. Since our model does not consider changes in density, it is not able to describe the effect of dilatancy and compaction, and thus under-pressure cannot be modeled, while the over-pressure phase, which is observed in our numerical simulations, is exclusively a consequence of advection (Fig. 7a-c). The relationship between distance along the x-axis and the maximum basal pore pressure reached is remarkably consistent with experimental data both in the curve shape and in the values measured (Fig. 7d), which further confirms the validity of the description of pore pressure used in our model once the under-pressure phase is finished, suggesting that the effect of compaction is limited compared to pore pressure advection. Note that the absence of dilatancy in our model is likely manifested in an earlier peak of the basal pore pressure than that expected in presence of an initial under-pressure stage (Fig. 7e).
Roche (2012) proposed that the basal under-pressure measured at the head of granular flows scales with the square of the flow front velocity. Based on this observation, Breard et al. (2019a) showed that the differential pressure measured in experiments can be given by

$$p_c = p_f - \frac{1}{2} \rho_b (u_{\text{slip}})^2, \qquad (11)$$

where $p_f$ is the basal pore pressure and $\rho_b$ is the mixture density at the base. The use of this expression and our numerical results show that the temporal evolution of $p_c$ at different positions along the channel is characterized by a short under-pressure phase followed by a longer over-pressure stage in proximal domains ($x < 0.5$ m, Fig. 8a-c), while distal points present only the under-pressure phase (Fig. 8a-c) because most of the pore pressure is already dissipated at these distances from the reservoir. Although the duration of the under-pressure phase of $p_c$ is shorter than that measured in the benchmark experiment, the simulated minimum values, their evolution with distance and the time at which these values are reached are strongly consistent with experimental data (Fig. 8d-e). On the other hand, while the maximum values of $p_c$ and those of experimental data are in agreement, the times at which these maximum values are observed are shifted.

The deposition dynamics of particles in the simulations is shown in Figure 9. Note that these results are a direct consequence of the rheological model adopted and no calibrated inputs of sedimentation rate are needed to parametrize the deposition of granular material. The length of the sliding head ($L_h$, Fig. 9a) was computed considering that sedimentation occurs at the base of the channel when the slip velocity reaches a value lower than 5% of the maximum slip velocity observed during the simulation. On the other hand, the variable $A_d$ (area of material deposited, Fig. 9a) was calculated by considering the modulus of velocity in each cell of the computational grid. In particular, at a given distance from the reservoir, the thickness of the deposit was computed considering all the cells with a velocity modulus lower than 0.1 m/s (i.e. about 5% of the maximum value of the velocity modulus observed during the simulations). Our
simulations show maximum lengths of the flow head of the order of 0.85-1.15 m, slightly larger than the experimental estimates of Roche (2012) (i.e. ~0.7 m; Fig. 9b) and twice the values simulated and observed in analogue experiments of dry granular flows of the same dimensions (i.e. 0.4 – 0.5 m; Roche 2012; Chupin et al. 2021). The relationship between $L_h/L$ and $L_d/L_f$ (see Fig. 9a for definition) shows a linear trend, in agreement with experimental data and also with the behaviour of dry granular flows (Roche 2012; Chupin et al. 2021; Fig. 9c). The evolution of the deposit area compared with the normalized time and run-out distance is also consistent with experimental data. Most of the deposition occurs during the final 40% of the propagation time-span, when the flow front has already travelled more than 80% of the final run-out distance (Fig. 9d-e). The results show that most of the deposition occurs between $t \approx 4.0\sqrt{H/g}$ and $t \approx 6.0\sqrt{H/g} - 6.5\sqrt{H/g}$ (Fig. 9f-g) and that lowering the pore pressure diffusion coefficient leads to delayed deposition. The position at which the peak of sedimentation rate occurs increases monotonically with time in all the simulations presented, showing slightly S-shaped curves that start in the vicinity of the reservoir and present maximum advance velocities similar in all the cases ($u_{sp}/\sqrt{gH} \sim 1.4$, where $u_{sp}$ is the advance velocity of the position of the deposition rate peak), significantly higher than the flow front velocity (Fig. 9f). Thus, our results suggest that the advance of the position of maximum deposition rate is poorly correlated with the behaviour of the flow front. At a given point along the channel, the deposition of particles tends to occur very rapidly (Fig. 9g). In fact, the time elapsed between the deposition of 10% and 90% of the final deposit at a given point is of the order of $0.1\sqrt{H/g}$, one order of magnitude smaller than the granular flow duration (Fig. 9g). Locally, the sedimentation rate reaches peaks of the order of 1 m/s, with mean sedimentation rates of the order of 0.1 m/s. It is worth noting, however, that this constraint of sedimentation rate is strongly influenced by the threshold used to define the deposited portion of the granular flow during its propagation.
385 **3.3 Simulations on inclined planes**

386 In the previous section we showed that the effective diffusion coefficient required to simulate
387 the benchmark analogue experiment (Roche et al. 2010), which is likely variable in time and
388 position, is in the range $\kappa = 0.015 - 0.035 \text{ m}^2/\text{s}$. Considering these values of diffusion
389 coefficient, we investigated the coupled effect of fluidization and topography through an
390 additional set of simulations adopting variable surface slope angles (from 0° to 10°).
391 Additionally, for comparison purposes, we did complementary dam-break simulations
392 considering inclined surfaces and dry flows, using the model described in Chupin et al. (2021).
393 Thereby, simulations results for run-out distance allow quantifying the combined effects of pore
394 pressure and surface slope angle.

395 The temporal evolution of the front position of dry and fluidized granular flows shows that a
396 small increment of the surface slope angle is able to significantly increase the maximum run-
397 out distance (Fig. 10). For instance, an increment of surface slope from 0° to 10° is able to
398 increase the modelled run-out distance from ~2H to ~4.4H for dry granular flows (relative
399 increase of 220%), where H is the initial column height, while an increase from ~3.8H to
400 ~8.4H was computed for fluidized flows with $\kappa = 0.025 \text{ m}^2/\text{s}$ (relative increase of 220%).
401 Differences in the propagation velocity between dry and fluidized granular flows are evident
402 from early phases of flow propagation (Figs. 10 and 11). On the other hand, as observed in the
403 simulations described in Section 3.2 for granular flows propagating on horizontal surfaces, also
404 on inclined planes the differences between simulations performed with different diffusion
405 coefficients are mainly manifested in the duration of the initial phase of velocity increase (and
406 thus manifested in the maximum front velocity reached by the flow; Fig. 11). Instead, for a
407 given slope angle, the velocity decrease during the final phase occurs at a similar rate for all the
408 fluidized and dry flows simulated (Fig. 11). This is because, once the pore pressure has been
409 dissipated by diffusion, the rheology of all the simulated granular flows is that of dry flows.
410

The deceleration of the flow front is strongly controlled by the slope angle (Fig. 11), ranging from ~0.23g (at 0°) to ~0.11g (at 10°). This dependency gives rise to significant differences in the modeled run-out distance as a function of surface slope angle for both dry and fluidized flows (Fig. 12).

411

412

413

414

Run-out of simulated flows also shows important aspects of pore pressure and surface slope angle effects (Fig. 12). For the range of diffusion coefficients adopted here, fluidization of the initial source is able to increase the run-out distance between ~1.55H ($\kappa$ = 0.035 m²/s, slope angle of 0°) and ~5.15H ($\kappa$ = 0.015 m²/s, slope angle of 10°), corresponding to an increase range of the run-out distance between ~165% and ~225%. Interestingly, the relative increase of run-out distance when fluidized granular flows are compared with dry flows is only weakly controlled by the slope angle (Fig. 12c). On the other hand, for given fluidization conditions, we note that an increase of slope angle from 0° to 10° produces an increment of the run-out distance of about 105-125%. This relative increase in the run-out distance is significantly larger than that measured in analogue experiments by Chédeville and Roche (2015) for lower-aspect ratio collapsing columns (0.5-1.0), i.e., ~60% for an increase of slope angle from 0° to 10°. We speculate that this could be a consequence of the slower pore pressure diffusion that characterizes taller collapsing columns.

415

416

417

418

419

420

421

422

423

424

425

426

427

The slope angle has a small influence on the maximum basal pore pressure computed at a given distance (Fig. 13). This shows that the evolution of the basal pore pressure is mainly controlled by the effective pore pressure diffusion coefficient. The length of the sliding head increases significantly when granular flows propagate on inclined surfaces (Fig. 14a-c). On the other hand, inclined topographies are able to delay the onset of deposition and reduce the sedimentation rate (Fig. 14d-f). Interestingly, the shape of the curves describing the evolution of the deposit area (Fig. 14d-f) changes when different slope angles are considered. Deposition in flows propagating on horizontal surfaces occurs at a nearly constant rate during almost all

428

429

430

431

432

433

434
the deposition stage (Fig. 14d-f), and the position at which the maximum deposition rate occurs
advances at an almost constant velocity (Fig. 14g-i). However, in simulations performed at high
slope angles, the initial stage of deposition, characterized by a relatively low sedimentation rate,
is accompanied by a relatively slow advance of the position at which the maximum deposition
rate occurs (Fig. 14g-i), while both the sedimentation rate and the advance velocity of the
position of maximum deposition increase during the final period of deposition (Fig. 14d-i).

441

4.42 **4. Discussion**

In this work we have presented a new model to describe dam-break fluidized granular flows
and test the effect of low-angle inclined surfaces in the resulting propagation dynamics. This
model, built on the formulation described by Chupin et al. (2021) for dry flows, was compared
with a benchmark analogue experiment for which detailed information of flow propagation,
pore pressure evolution and sedimentation dynamics is available in the literature (Roche et al.
2010; Roche 2012). Thereby, this work complements previous efforts to analyse analogue
experiments through numerical modeling (Breard et al. 2019a). In particular, Breard et al.
(2019a) tested different friction models and compared their simulations with experiments
considering flow shape, kinematics and pore pressure evolution. Our model allows the
description of the sedimentation dynamics of granular flows and their comparison with
additional characteristics of the benchmark experiment (Roche et al. 2010; Roche 2012), thus
allowing to explore aspects of granular flows that were not addressed by Breard et al. (2019a).
In contrast, the model of Breard et al. (2019a) is able to describe slight compaction and dilation
processes, which is not possible in our formulation.

Numerical results reproduce reasonably well the collapse and propagation dynamics described
in the analogue experiment in terms of run-out distance and pore pressure, and they allow to
constrain the effective diffusion coefficient that characterizes the granular material considered.
However, even though the model captures the general shape of the resulting deposits, the thickness tends to be under-estimated in proximal domains and over-estimated in distal domains. Potential sources of systematic differences between analogue experiments and our numerical model are the dynamics of gate opening and simplifications in the mathematical description such as the non-compressibility of the granular flow and the assumption of a constant effective diffusion coefficient in space and time.

Interestingly, our estimates of the effective diffusion coefficient are consistent with experimental measurements on static defluidizing beds (Roche et al. 2010) and are one order of magnitude smaller than the theoretical value, which fails completely in predicting the behaviour of the studied analogue experiment (see Fig. 1d). The discrepancy between the theoretical value and experiment-derived estimates (Roche et al. 2010; Montserrat et al. 2012) is a major unsolved issue related to pore pressure diffusion in granular materials. Breard et al. (2019b) showed that if the volume of air in a windbox at the base of an experimental granular column is significant compared to the volume of air in the column, then the measured diffusion coefficient is larger than predicted theoretically. However, we made recently further pore pressure diffusion tests in a device with a windbox whose volume was less than ~0.05 % of the volume of air in the granular column, and we found a positive correlation between the diffusion coefficient and the column height (in preparation). Therefore, though a windbox affects the estimates of pore pressure diffusion coefficient, it cannot be invoked to explain differences of more than one order of magnitude between experimental and theoretical estimates, and thus additional investigation is still required to understand this discrepancy. In the case of the numerical simulations presented here, we note that the effective diffusion coefficients $\kappa = 0.015 - 0.035 \text{ m}^2/\text{s}$ giving the best agreement with the experimental data are those typical of static bed heights of ~15–25 cm, which are about half the height of the initial column in the
484

dam-break configuration. This typical height seems to be the best compromise between the
height of the column released and that of the resulting flow.

485

486

Despite that main differences in flow dynamics due to different diffusion coefficients arise
during only the first ~17-22% of the total propagation time, they can cause significant changes
in the resulting run-out distance. In contrast, during the later phases of flow propagation, once
pore pressure has diffused significantly, the non-fluidized conditions of the flow produce a
similar stopping dynamics in all the simulations studied. These results suggest that
understanding the processes controlling the generation and evolution of pore pressure (e.g.
internal gas-particle motion, air ingestion, particle settling and diffusion; Sweeney and
Valentine 2017; Valentine and Sweeney 2018; Valentine 2020; Fries et al. 2021) at early
propagation stages can be particularly critical in controlling the whole granular flow, regardless
of possible mechanisms able to generate pore pressure during later propagation stages (Benage
et al. 2016; Breard et al. 2018; Chédeville and Roche 2015, 2018; Lube et al. 2019), which are
not taken into account in our numerical model and are expected to be negligible in the
benchmark experiment. In simulations on horizontal surfaces with effective diffusion
coefficients compatible with the benchmark experiment, we observe an increase of run-out
distance by a factor of ~1.8-2.2 when compared with dry granular flows. Thus, fluidization
processes represent a critical factor in the evaluation of PDC hazard.

487

488

489

490

491

492

493

494

495

496

497

498

499

500

501

502

Additionally, this work provides insights for understanding some aspects of the dynamics of
fluidized granular flows such as the evolution of pore pressure in time and space, the deposition
process, and the effect of inclined topographies. These aspects are discussed below:

503

504

505

a) Our simulations of initially fluidized flows present an initial phase dominated by pore pressure advection and a later phase controlled by pore pressure diffusion up to reach stationary conditions. The transition between these phases is influenced by the effect of front velocity on flow stretching because pore pressure diffuses faster in thinner flows.

506

507

508
509

510

511

512

513

Importantly, these results suggest that the fluidization effect in increasing the maximum run-out distance may be self-limited, particularly on steep slopes. In fact, high pore pressure reduces friction and causes faster granular flows able to travel larger distances, but in turn fast propagation causes reduction in flow thickness, which causes faster pore pressure diffusion.

514

515

516

517

518

519

520

521

522

523

524

525

526

527

528

529

530

531

b) The basal pore pressure simulated at a given point along the channel shows an over-pressure phase coincident with the passage of the flow head, while the pressure signal measured in experiments beneath the flow substrate interface, which represents the difference between the pressure at the flow base and the ambient atmospheric pressure, is characterized by a short under-pressure phase followed by a longer over-pressure stage. Although comparing these data is not straightforward because the experimental data are partially influenced by processes not considered in our model, in this work we analyzed the main characteristics of these signals and we also considered the influence of the basal slip conditions, as suggested in the literature (Roche et al. 2010; Breard et al. 2019a). We show that the relationship between distance along the channel and the maximum pressure reached during the flow passage is remarkably similar in simulations and experiments, which indicates that our model is able to capture reasonably well the evolution of pore pressure within the granular flow. This suggests that the effect of compaction and dilatancy processes (Bouchut et al. 2016, 2021) is limited once the flow front has passed, and that the pore pressure effect in the propagation of granular flows can be modeled considering only advection. Moreover, we show that the magnitude of the under-pressure phase measured in experiments can be successfully quantified by considering the slip velocity at the channel base, as proposed by Breard et al. (2019a).

532

533

c) Our simulations suggest that deposition is close to the *en masse* end-member. In fact, for a given point along the channel, the time-span during which deposition occurs is
534

535

536

537

538

539

540

541

much smaller than the timescale of granular flow propagation. Our results show that the
position at which the maximum of sedimentation rate occurs advances monotonically
from the reservoir and it is strongly influenced by the surface slope angle, while the
effect of the pore pressure diffusion coefficient is small. Our conclusion on deposition
is relevant for the experimental configuration considered but it is not necessarily
applicable for natural systems of significantly larger scale. In nature, in fact, progressive
aggradation can operate if onset of deposition occurs at early stages, the flow thickness
is large and the deposition rate is low (see Fig. 12 of Roche 2012).

542

543

544

545

546

547

548

549

550

551

552

553

554

555

556

557

558

d) Numerical simulations on inclined surfaces have shown that a low slope angle (up to 10°) is able to increase the run-out distance by a factor of 2.05 – 2.25 when compared with horizontal surfaces. This has major implications for pyroclastic density currents, which typically propagate at gentle slope angles. A remarkable example where the regional slope could exert a significant effect is the Cerro Galan Ignimbrite (NW Argentina; Francis et al. 1983; Cas et al. 2011; Lesti et al. 2011), which presents a maximum run-out distance of ~70 km and was emplaced on a regional regular slope of a few degrees. Cas et al. (2011) also suggested an important effect of gas pore pressure in the reduction of friction between the flow and the substrate in this case study. Another example is the Peach Spring Tuff (USA), formed by PDCs that travelled >170 km from the eruptive centre and propagated on substrates with gentle slope angles (Valentine et al. 1989; Roche et al. 2016). In this case as well, a significant influence of gas pore pressure in the resulting run-out distance has been suggested (Roche et al. 2016). Notice that though regional slope may enhance the runout distance of PDCs, recent advances suggest that the latter is controlled fundamentally by the discharge rate (Roche et al. 2021).
559 **5. Concluding remarks**
560

The numerical simulations presented in this work and their comparison with published
experimental data have revealed that:

562 (1) Even though the pore pressure diffusion coefficient probably varies in space and time
563 in dam-break fluidized granular flows, a constant (effective) pore pressure diffusion
564 coefficient can be estimated to capture reasonably well the flow dynamics in terms of
565 run-out distance, temporal evolution of pore fluid pressure, and shape of the deposit.
566

566 (2) Pore pressure increases significantly the run-out distance of initially fluidized granular
567 flows when compared with dry granular flows (e.g., by a factor of ~1.8 – ~2.2 on
568 horizontal slopes). Therefore, taking into account pore fluid pressure appears critical for
569 modelling dense PDCs in the context of volcanic hazard assessment.
570

570 (3) A significant effect in granular flow run-out is also exerted by the substrate slope angle.
571 For given fluidization conditions, an increase of slope angle from 0° to 10° produces an
572 increment of the run-out distance of 105 – 125%.
573

573 (4) The effect of fluidization in increasing run-out distance may be self-limited because the
574 higher velocity due to fluidization tends to reduce flow thickness, which induces faster
575 pore pressure diffusion.
576

576 (5) The pore pressure evolution in initially fluidized granular flows is mainly controlled by
577 the diffusion coefficient, while the effect of the angle slope of the substrate is limited.
578

578 (6) In the dam-break configuration at laboratory scale, the onset of the deposition of
579 granular flows occurs with a significant delay with respect to the front propagation.
580 Once deposition starts, the position at which the maximum sedimentation rate occurs
581 advances monotonically with time at a velocity significantly larger than the flow front
582 velocity. The dynamics of sedimentation in the studied experimental configuration,
583 which is a direct consequence of the rheological model adopted and does not require
584

calibrated inputs to set the sedimentation rate, is close to the *en masse* end-member model, but more progressive aggradation may operate in nature.

585

586

(7) Our model describes depth-dependent variations of the properties of granular flows considering high-aspect ratio dam-break configurations. Moreover, the possibility of exploring granular flows at larger length-scale makes this model a promising tool for investigating the factors controlling the dynamics of long run-out PDCs in nature.

587

588

589

590

591

## Acknowledgement

592

This research was financed by the French government IDEX-ISITE initiative 16-IDEX-0001 (CAP 20-25). This is Laboratory of Excellence ClerVolc contribution number 501. The numerical simulations have been performed on a DELL cluster with 32 processors Xeon E2650v2 (8 cores), 1 To of total memory and an infiniband (FDR 56Gb/s) connecting network.

593

594

595

596

This cluster has been financed by the French Government Laboratory of Excellence initiative n°ANR-10-LABX-0006. We thank an anonymous referee, Dr. Gert Lube and Dr. Greg Valentine for their useful comments and suggestions to improve this work.

597

598

599

600

## References

601

Balay S, Abhyankar S, Adams M, Brown J, Brune P, Buschelman K, Dalcin L, Dener A, Eijkhout V, Gropp W, Karpeyev D, Kaushik D, Knepley M, May D, McInnes LC, Mills R, Munson T, Rupp K, Sanan P, Smith B, Zampini S, Zhang H (2018a) PETSc Users Manual Revision 3.10 (Technical Report ANL-95/11 - Revision 3.10). Argonne National Laboratory (ANL), Argonne, IL (United States).

602

603

604

605

606

Balay S, Abhyankar S, Adams M, Brown J, Brune P, Buschelman K, Dalcin L, Dener A, Eijkhout V, Gropp W, Kaushik D, Knepley M, May D, McInnes LC, Mills R, Munson

607
608

T, Rupp K, Sanan P, Smith B, Zampini S, Zhang H (2018b) PETSc Web page.

609

http://www.mcs.anl.gov/petsc

610

Benage M, Dufek J, Mothes P (2016) Quantifying entrainment in pyroclastic density currents

611

from the Tungurahua eruption, Ecuador: Integrating field proxies with numerical

612

simulations. Geophysical Research Letters 43(13):6932-6941.

613

https://doi.org/10.1002/2016GL069527

614

Bernard J, Kelfoun K, Le Pennec J, Vallejo Vargas S (2014) Pyroclastic flow erosion and

615

bulking processes: comparing field-based vs modeling results at Tungurahua volcano,

616

Ecuador. Bulletin of Volcanology 76(9):1-16.

617

https://doi.org/10.1007/s00445-014-0858-y

618

Bouchut F, Fernández-Nieto E, Mangeney A, Narbona-Reina G (2016) A two-phase two-layer

619

model for fluidized granular flows with dilatancy effects. Journal of Fluid Mechanics

620

801:166-221. https://doi.org/10.1017/jfm.2016.417

621

Bouchut F, Fernández-Nieto E, Mangeney A, Narbona-Reina G (2021) Dilatancy in dry

622

granular flows with a compressible µ(I) rheology. Journal of Computational Physics

623

429:110013. https://doi.org/10.1016/j.jcp.2020.110013

624

Brand B, Mackaman-Lofland C, Pollock N, Bendaña S, Dawson B, Wichgers P (2014)

625

Dynamics of pyroclastic density currents: Conditions that promote substrate erosion and

626

self-channelization—Mount St Helens, Washington (USA). Journal of Volcanology and

627

Geothermal Research 276:189-214. https://doi.org/10.1016/j.jvolgeores.2014.01.007

628

Branney M, Kokelaar P (2002) Pyroclastic density currents and the sedimentation of

629

ignimbrites. Geological Society of London.

630

Breard EC, Lube G, Jones JR, Dufek J, Cronin SJ, Valentine GA, Moebis A (2016) Coupling

631

of turbulent and non-turbulent flow regimes within pyroclastic density currents. Nature

632

Geoscience 9(10):767-771. https://doi.org/10.1038/ngeo2794
633

Breard E, Dufek J, Lube G (2018) Enhanced mobility in concentrated pyroclastic density currents: An examination of a self-fluidization mechanism. Geophysical Research Letters, 45(2), 654-664. https://doi.org/10.1002/2017GL075759

634

635

636

Breard EC, Dufek J, Roche O (2019a) Continuum modeling of pressure-balanced and fluidized granular flows in 2-D: Comparison with glass bead experiments and implications for concentrated pyroclastic density currents. Journal of Geophysical Research: Solid Earth 124(6):5557-5583. https://doi.org/10.1029/2018JB016874

637

638

639

640

Breard EC, Jones, JR, Fullard L, Lube G, Davies C, Dufek J (2019b) The permeability of volcanic mixtures—implications for pyroclastic currents. Journal of Geophysical Research: Solid Earth 124(2):1343-1360. https://doi.org/10.1029/2018JB016544

641

642

643

Burgisser A (2012) A semi-empirical method to calculate the permeability of homogeneously fluidized pyroclastic material. Journal of Volcanology and Geothermal Research 243: 97-106. https://doi.org/10.1016/j.jvolgeores.2012.08.015

644

645

646

Burgisser A, Bergantz G (2002) Reconciling pyroclastic flow and surge: the multiphase physics of pyroclastic density currents. Earth and Planetary Science Letters 202(2):405-418. https://doi.org/10.1016/S0012-821X(02)00789-6

647

648

649

Bursik MI, Woods AW (1996) The dynamics and thermodynamics of large ash-flows. Bulletin of Volcanology 58:175-93. https://doi.org/10.1007/s004450050134

650

651

Cas R, Wright H, Folkes C, Lesti C, Porreca M, Giordano G, Viramonte J (2011) The flow dynamics of an extremely large volume pyroclastic flow, the 2.08-Ma Cerro Galán Ignimbrite, NW Argentina, and comparison with other flow types. Bulletin of Volcanology 73(10):1583-1609. https://doi.org/10.1007/s00445-011-0564-y

652

653

654

655

Chalayer R, Chupin L, Dubois T (2018) A bi-projection method for incompressible bingham flows with variable density, viscosity, and yield stress. SIAM Journal on Numerical Analysis 56(4):2461-2483. https://doi.org/10.1137/17M113993X

656

657
658 Chédeville C, Roche O (2015) Influence of slope angle on pore pressure generation and kinematics of pyroclastic flows: insights from laboratory experiments. Bulletin of Volcanology 77(11):1-13. https://doi.org/10.1007/s00445-015-0981-4

661 Chédeville C, Roche O (2018) Autofluidization of collapsing bed of fine particles: implications for the emplacement of pyroclastic flows. Journal of Volcanology and Geothermal Research 368:91-99. https://doi.org/10.1016/j.jvolgeores.2018.11.007

664 Chupin L, Dubois T (2016) A bi-projection method for Bingham type flows. Computers & Mathematics with Applications 72(5):1263-1286. https://doi.org/10.1016/j.camwa.2016.06.026

667 Chupin L, Dubois T, Phan M, Roche O (2021) Pressure-dependent threshold in a granular flow: Numerical modeling and experimental validation. Journal of Non-Newtonian Fluid Mechanics 291:104529. https://doi.org/10.1016/j.jnnfm.2021.104529

670 Cole P, Neri A, Baxter P (2015) Hazards from pyroclastic density currents. In: The Encyclopedia of Volcanoes (pp. 943-956). Academic Press.

672 Druitt T (1998) Pyroclastic density currents. Geological Society, London, Special Publications 145(1):145-182.

674 Druitt T, Avard G, Bruni G, Lettieri P, Maez F (2007) Gas retention in fine-grained pyroclastic flow materials at high temperatures. Bulletin of Volcanology 69(8):881-901. https://doi.org/10.1007/s00445-007-0116-7

677 Dufek J (2016) The fluid mechanics of pyroclastic density currents. Annual Review of Fluid Mechanics 48:459-485. https://doi.org/10.1146/annurev-fluid-122414-034252

679 Dufek J, Esposti Ongaro T, Roche O (2015) Pyroclastic density currents: processes and models. In: The Encyclopedia of Volcanoes (pp. 617-629). Academic Press.
681

Esposti Ongaro T, Cerminara M, Charbonnier S, Lube G, Valentine G (2020) A framework for validation and benchmarking of pyroclastic current models. Bulletin of Volcanology 82:1-17. https://doi.org/10.1007/s00445-020-01388-2

682

683

684

Farin M, Mangeney A, Roche O (2014) Fundamental changes of granular flow dynamics, deposition, and erosion processes at high slope angles: insights from laboratory experiments. Journal of Geophysical Research: Earth Surface 119(3):504-532. https://doi.org/10.1002/2013JF002750

685

686

687

688

Francis P, O'Callaghan L, Kretzschmar G, Thorpe R, Sparks R, Page R, De Barrio R, Gillou G, Gonzalez O (1983) The Cerro Galan ignimbrite. Nature 301(5895):51-53. https://doi.org/10.1038/301051a0

689

690

691

Fries A, Roche O, Carazzo G (2021) Granular mixture deflation and generation of pore fluid pressure at the impact zone of a pyroclastic fountain: Experimental insights. Journal of Volcanology and Geothermal Research 414:107226. https://doi.org/10.1016/j.jvolgeores.2021.107226

692

693

694

695

Fullmer WD, Hrenya CM (2017) The clustering instability in rapid granular and gas-solid flows. Annual Review of Fluid Mechanics 49:485-510. https://doi.org/10.1146/annurev-fluid-010816-060028

696

697

698

Garres-Díaz J, Bouchut F, Fernández-Nieto E, Mangeney A, Narbona-Reina G (2020) Multilayer models for shallow two-phase debris flows with dilatancy effects. Journal of Computational Physics 419:109699. https://doi.org/10.1016/j.jcp.2020.109699

699

700

701

Gase A, Brand B, Bradford J (2017) Evidence of erosional self-channelization of pyroclastic density currents revealed by ground-penetrating radar imaging at Mount St. Helens, Washington (USA). Geophysical Research Letters 44(5):2220-2228. https://doi.org/10.1002/2016GL072178

702

703

704
705

Giordano G, Cas RAF (2021) Classification of ignimbrites and their eruptions. Earth-Science

706

Reviews 220:103697. https://doi.org/10.1016/j.earscirev.2021.103697

707

Goren L, Aharonov E, Sparks D, Toussaint R (2010) Pore pressure evolution in deforming

708

granular material: A general formulation and the infinitely stiff approximation. Journal

709

of Geophysical Research: Solid Earth 115(B9). https://doi.org/10.1029/2009JB007191

710

Gray J, Edwards A (2014) A depth-averaged µ(I)-rheology for shallow granular free-surface

711

flows. Journal of Fluid Mechanics 755:503-534. https://doi.org/10.1017/jfm.2014.450

712

Gueugneau V, Kelfoun K, Roche O, Chupin L (2017) Effects of pore pressure in pyroclastic

713

flows: numerical simulation and experimental validation. Geophysical Research Letters

714

44(5):2194-2202. https://doi.org/10.1002/2017GL072591

715

Ionescu I, Mangeney A, Bouchut F, Roche O (2015) Viscoplastic modeling of granular column

716

collapse with pressure-dependent rheology. Journal of Non-Newtonian Fluid Mechanics

717

219:1-18. https://doi.org/10.1016/j.jnnfm.2015.02.006

718

Iverson R (1997) The physics of debris flows. Reviews of Geophysics 35(3):245-296.

719

https://doi.org/10.1029/97RG00426

720

Iverson R, Denlinger R (2001) Flow of variably fluidized granular masses across three-

721

dimensional terrain: 1. Coulomb mixture theory. Journal of Geophysical Research:

722

Solid Earth 106(B1):537-552. https://doi.org/10.1029/2000JB900329

723

Jop P, Forterre Y, Pouliquen O (2006) A constitutive law for dense granular flows. Nature

724

441(7094):727-730. https://doi.org/10.1038/nature04801

725

Kelfoun K (2011) Suitability of simple rheological laws for the numerical simulation of dense

726

pyroclastic flows and long-runout volcanic avalanches. Journal of Geophysical

727

Research: Solid Earth 116(B8). https://doi.org/10.1029/2010JB007622

728

Lesti C, Porreca M, Giordano G, Mattei M, Cas R, Wright H, Folkes C, Viramonte J (2011)

729

High-temperature emplacement of the Cerro Galán and Toconquis Group ignimbrites
730 (Puna plateau, NW Argentina) determined by TRM analyses. Bulletin of Volcanology
731 73(10):1535-1565. https://doi.org/10.1007/s00445-011-0536-2

732 Lube G, Breard EC, Jones J, Fullard L, Dufek J, Cronin SJ, Wang T (2019) Generation of air
733 lubrication within pyroclastic density currents. Nature Geoscience 12(5):381-386.
734 https://doi.org/10.1038/s41561-019-0338-2

735 Lube G, Breard EC, Esposti Ongaro T, Dufek J, Brand B (2020) Multiphase flow behaviour
736 and hazard prediction of pyroclastic density currents. Nature Reviews Earth &
737 Environment 1(7):348-365. https://doi.org/10.1038/s43017-020-0064-8

738 Marino B, Thomas L, Linden P (2005) The front condition for gravity currents. Journal of Fluid
739 Mechanics 536:49-78. https://doi.org/10.1017/S0022112005004933

740 Montserrat S, Tamburrino A, Roche O, Niño Y (2012) Pore fluid pressure diffusion in
defluidizing granular columns. Journal of Geophysical Research: Earth Surface
741 117(F2). https://doi.org/10.1029/2011JF002164

743 Neri A, Esposti Ongaro T, Voight B, Widiwijayanti C (2015) Pyroclastic density current
744 hazards and risk. In: Volcanic Hazards, Risks and Disasters (pp. 109-140). Elsevier.

745 Osher S, Fedkiw R (2001) Level set methods: an overview and some recent results. Journal of
746 Computational Physics 169(2):463-502. https://doi.org/10.1006/jcph.2000.6636

747 Roche O (2012) Depositional processes and gas pore pressure in pyroclastic flows: an
748 experimental perspective. Bulletin of Volcanology 74(8):1807-1820.
749 https://doi.org/10.1007/s00445-012-0639-4

750 Roche O, Montserrat S, Niño Y, Tamburrino A (2010) Pore fluid pressure and internal
751 kinematics of gravitational laboratory air-particle flows: Insights into the emplacement
752 dynamics of pyroclastic flows. Journal of Geophysical Research: Solid Earth 115(B9).
753 https://doi.org/10.1029/2009JB007133
754

Roche O, Buesch D, Valentine G (2016) Slow-moving and far-travelled dense pyroclastic flows

755

during the Peach Spring super-eruption. Nature Communications 7(1):1-8.

756

https://doi.org/10.1038/ncomms10890

757

Roche O, Azzaoui N, Guillin A (2021) Discharge rate of explosive volcanic eruption controls

758

runout distance of pyroclastic density currents. Earth and Planetary Science Letters

759

568:117017. https://doi.org/10.1016/j.epsl.2021.117017

760

Rowley P, Roche O, Druitt T, Cas R (2014) Experimental study of dense pyroclastic density

761

currents using sustained, gas-fluidized granular flows. Bulletin of Volcanology 76(9):1-

762

13. https://doi.org/10.1007/s00445-014-0855-1

763

Savage S, Iverson R (2003) Surge dynamics coupled to pore-pressure evolution in debris flows.

764

In: Proc. 3rd Int. Conf. on Debris-Flow Hazards Mitigation: Mechanics, Prediction, and

765

Assessment. Edited by: Rickenmann D, Chen CL. Millpress, Rotterdam (pp. 503-514).

766

Shimizu H, Koyaguchi T, Suzuki Y (2019) The run-out distance of large-scale pyroclastic

767

density currents: a two-layer depth-averaged model. Journal of Volcanology and

768

Geothermal Research 381:168-184. https://doi.org/10.1016/j.jvolgeores.2019.03.013

769

Smagorinsky J (1963) General circulation experiments with the primitive equations: I. The

770

basic experiment. Monthly Weather Review 91(3):99-164.

771

Sweeney M, Valentine G (2017) Impact zone dynamics of dilute mono-and polydisperse jets

772

and their implications for the initial conditions of pyroclastic density currents. Physics

773

of Fluids 29(9):093304. https://doi.org/10.1063/1.5004197

774

Valentine G (2019) Preface to the topical collection—pyroclastic current models:

775

benchmarking and validation. Bulletin of Volcanology 81:69.

776

https://doi.org/10.1007/s00445-019-1328-3
777

Valentine G (2020) Initiation of dilute and concentrated pyroclastic currents from collapsing

778

mixtures and origin of their proximal deposits. Bulletin of Volcanology 82(2):1-24.

779

https://doi.org/10.1007/s00445-020-1366-x

780

Valentine G, Buesch D, Fisher R (1989) Basal layered deposits of the Peach Springs Tuff,

781

northwestern Arizona, USA. Bulletin of Volcanology 51(6):395-414.

782

https://doi.org/10.1007/BF01078808

783

Valentine G, Doronzo D, Dellino P, de Tullio M (2011) Effects of volcano profile on dilute

784

pyroclastic density currents: numerical simulations. Geology 39(10):947-950.

785

https://doi.org/10.1130/G31936.1

786

Valentine G, Sweeney M (2018) Compressible flow phenomena at inception of lateral density

787

currents fed by collapsing gas-particle mixtures. Journal of Geophysical Research: Solid

788

Earth 123(2):1286-1302. https://doi.org/10.1002/2017JB015129

789

Wilson C (1984) The role of fluidization in the emplacement of pyroclastic flows, 2:

790

experimental results and their interpretation. Journal of Volcanology and Geothermal

791

Research 20(1-2):55-84. https://doi.org/10.1016/0377-0273(84)90066-0

792
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Surface profiles of the granular flows at different times after release (see legends) in four simulations performed on horizontal planes, considering initially fluidized conditions and different values of the effective diffusion coefficient ($\kappa$, see titles). The final surface profile of the benchmark analogue experiment is also included (Roche et al. 2010).</figcaption></figure>

794

795

796

797

798

799
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Surface profiles of the granular flows at different times after release in four simulations performed on horizontal planes, considering initially fluidized conditions and different values of the effective diffusion coefficient ($\kappa$, see legend). The evolution of the surface profile of the benchmark analogue experiment is also included (Roche et al. 2010).</figcaption></figure>

800

801

802

803

804

805

806
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Temporal evolution of the front position of the granular flows in four simulations performed on horizontal planes, considering initially fluidized conditions and different values of the effective diffusion coefficient ($κ$, see legend). The temporal evolution of the front position of the benchmark analogue experiment is shown as well (see legend; Roche et al. 2010). Because our simulations are not able to describe flow thicknesses lower than 3.5 mm (i.e. the cell size used in numerical simulations), we also include the experimental data considering only flow thicknesses above this threshold in the definition of the front position. Both axes are normalized using ad-hoc factors in order to produce non-dimensional results ($H = 0.4$ m and $g = 9.8$ m/s$^2$).</figcaption></figure>

807

808

809

810

811

812

813

814

815

816

817

818
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Temporal evolution of the front velocity of the granular flows in four simulations performed on horizontal planes, considering initially fluidized conditions and different values of the effective diffusion coefficient ($κ$, see legend). A moving average function was applied to these curves, considering a time window of 0.1 s. The evolution of the front velocity of the benchmark analogue experiment is shown as well (see legend; Roche et al. 2010). Because our simulations are not able to describe flow thicknesses lower than 3.5 mm (i.e. the cell size used in numerical simulations), we include the experimental data considering only flow thicknesses above this threshold in the definition of the front position. Both axes are normalized using ad-hoc factors in order to produce non-dimensional results ($H = 0.4$ m and $g = 9.8$ m/s²). Note that the theoretical value for the maximum velocity in a dam-break experiment of an inviscid flow is $u/\sqrt{gH} = \sqrt{2} \approx 1.4$ (Marino et al. 2005).</figcaption></figure>

819

820

821

822

823

824

825

826

827

828

829

830

831

832
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Basal pore pressure profiles of the granular flows at different times after release (see legends) in three simulations performed on horizontal planes, considering initially fluidized conditions and different values of the effective diffusion coefficient ($κ$, see titles). Note that the ratio between basal pore pressure and the lithostatic pressure ($p_f/ρgh$), not shown here, represents the degree of fluidization (see Supplementary Figure S1). Full fluidization occurs when $p_f/ρgh$ is larger than one.</figcaption></figure>

833

834

835

836

837

838

839

840

841
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Temporal evolution of the front position and a set of iso-pressure fronts (i.e. the position along the x-axis at which specific values of pore pressure are reached as a function of time, see legends) in three simulations performed on horizontal planes, considering initially fluidized conditions and different values of the effective diffusion coefficient ($κ$, see titles).</figcaption></figure>

842

843

844

845

846

847

848
<figure><img src="image_9.png" /><figcaption>**Figure 7.** (a-c) Temporal evolution of basal pore pressure at different positions (see legends) in simulations on horizontal planes, considering initially fluidized conditions and different values of the effective diffusion coefficient ($\kappa$, see titles). Experimental data are also presented (Roche et al. 2010). (d) Maximum normalized values of basal pressure in numerical simulations and in the benchmark experiment as a function of horizontal distance (see legend; Roche et al. 2010). (e) Time needed to reach the extreme values of basal pressure (i.e. minimum, if present, and maximum values) at the channel base in numerical simulations and in the benchmark experiment (see legend; Roche et al. 2010).</figcaption></figure>

849

850

851

852

853

854

855

856

857
<figure><img src="image_10.png" /><figcaption>**Figure 8.** (a-c) Temporal evolution of $p_c$ (see equation (11)) at different positions (see legends) in simulations performed on horizontal planes, considering initially fluidized conditions and different values of the effective diffusion coefficient ($\kappa$, see titles). Experimental data are also presented (Roche et al. 2010), which describe the difference between the pressure generated by the flow above a sensor located at $x = 0.2$ m and the ambient atmospheric pressure. (d) Maximum normalized values of $p_c$ and differential pressure with respect to the atmosphere in numerical simulations and in the benchmark experiment, respectively, as a function of</figcaption></figure>

858

859

860

861

862

863

864

865
866

horizontal distance (see legend; Roche et al. 2010). (e) Time needed to reach the extreme values

867

(i.e. minimum and maximum values) of $p_c$ and of differential pressure with respect to the

868

atmosphere in numerical simulations and in the benchmark experiment (Roche et al. 2010),

869

respectively (see legend).

870
<figure><img src="image_11.png" /></figure>

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

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

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

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

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

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

871

872

<caption>**Figure 9.** Plots describing the deposition dynamics of particles in our numerical simulations.</caption>

873

874

875

(a) Schematic figure showing the definitions used to describe the deposition dynamics of the modeled granular flows. (b) Temporal evolution of the sliding head length ($L_h$, see panel a) in simulations performed on horizontal planes, considering initially fluidized conditions and
876

different values of the effective diffusion coefficient ($\kappa$, see legend). (c) $L_h/L$ as a function of $L_d/L_f$ (see panel a) in the same set of simulations, where $L_f = L(t_f)$ and $t_f$ is the final time.

877

878

(d) Temporal evolution of $A_d/A_d(t_f)$ (see panel a) in the same set of simulations. (e) $A_d/A_d(t_f)$ as a function of $L/L_f$ (see panel a) in the same set of simulations. (f) Temporal

879

880

881

882

883

884

evolution of the position at which the peak of deposition rate is modeled in the same set of simulations. The front position is also included. (g) Times at which different percentiles (10%, 50% and 90%) of the final deposit thickness are reached as a function of distance along the x-axis, for the same set of simulations. In panels (b)-(e) we include data from the benchmark analogue experiment (Roche 2012).

885

886
<figure><img src="image_18.png" /><figcaption>(a)</figcaption></figure>

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

<figure><img src="image_20.png" /><figcaption>(c)</figcaption></figure>

887

888

889

890

891

892

<caption>**Figure 10.** Temporal evolution of the front position of the granular flows in simulations with variable initial fluidization conditions (dry and fluidized flows) and different values of the effective diffusion coefficient and surface slope angle (see titles and legend).</caption>
<figure><img src="image_21.png" /><figcaption>**Figure 11.** Temporal evolution of the front velocity of the granular flows in simulations with variable initial fluidization conditions (dry and fluidized flows) and different values of the effective diffusion coefficient and surface slope angle (see titles and legend). A moving average function was applied to these curves, considering a time window of 0.1 s.</figcaption></figure>

893

894

895

896

897

898

899
<figure><img src="image_22.png" /><figcaption>**Figure 12.** Comparison of the run-out distance of granular flows in simulations performed considering variable initial conditions (dry and fluidized flows) and different effective diffusion coefficients as a function of the surface slope angle. We use the function $\sin(\cdot)$ in the x-axes because it is the driving component of gravity. RD/H: run-out distance over initial column height (H). (RD – RD<sub>dry</sub>)/H: increase of run-out distance over H with respect to dry flows. RD/RD<sub>dry</sub>: ratio of run-out distance with respect to dry flows. (RD – RD<sub>hor</sub>)/H: increase of run-out distance over H with respect to a flow propagated over a horizontal surface. RD/RD<sub>hor</sub>: ratio of run-out distance with respect to a flow propagated over a horizontal surface.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Figure 13.** Maximum normalized basal pore pressure during the propagation of fluidized granular flows as a function of horizontal distance in numerical simulations with different effective diffusion coefficients and surface slope angles (see legend).</figcaption></figure>

911

912

913

914

915

916
<figure><img src="image_24.png" /><figcaption>**Figure 14.** (a-c) Temporal evolution of $L_h$ (see Figure 9a) in simulations performed on horizontal and inclined planes (see legends), considering initially fluidized conditions and different values of the effective diffusion coefficient ($\kappa$, see titles). (d-f) Temporal evolution of $A_d/A_d(t_f)$ (see Figure 9a) in the same set of simulations. (g-i) Temporal evolution of the position at which the peak of deposition rate is modeled in the same set of simulations. The front position is also included.</figcaption></figure>

917

918

919

920

921

922

923

924