# Modelling Settling-Driven Gravitational Instabilities at the Base of
Volcanic Clouds Using the Lattice Boltzmann Method

Jonathan Lemus1,2*, Allan Fries1, Paul A. Jarvis1, Costanza Bonadonna1,
Bastien Chopard2, Jonas Lätt2
1Department of Earth Sciences, University of Geneva, 1205 Geneva, Switzerland
2Department of Computer Science, University of Geneva, 1227 Carouge,
Switzerland
<EMAIL_ADDRESS>
## 1\. Introduction

Explosive volcanic eruptions can inject large quantities of ash into the
atmosphere, generating multiple hazards at various spatial and temporal scales
(Blong,, 2000; Bonadonna et al.,, 2021). Subsequent volcanic ash dispersal and
sedimentation can strongly disrupt air traffic (Guffanti et al.,, 2008; Prata
and Rose,, 2015), affect inhabited areas (Jenkins et al.,, 2015; Spence et
al.,, 2005), and impact ecosystems and public health (Gudmundsson,, 2011;
Wilson et al.,, 2011). A good understanding of ash dispersal is critical for
effective forecasting and management of the response to these hazards. Modern
volcanic ash transport and dispersal models have now reached high levels of
sophistication (Bonadonna et al.,, 2012; Folch,, 2012; Folch et al.,, 2020;
Jones et al.,, 2007; Prata et al.,, 2021) but do not include all of the
physical processes affecting ash transport, such as particle aggregation and
settling-driven gravitational instabilities (e.g. Durant,, 2015). Various
studies have highlighted the need to take these processes into account by
revealing discrepancies between field measurements and numerical models
(Scollo et al.,, 2008), premature sedimentation of fine ash leading to bimodal
grainsize distributions not only related to particle aggregation (Bonadonna et
al.,, 2011; Manzella et al.,, 2015; Watt et al.,, 2015) and significant
depletion of airborne fine ash close to the source (Gouhier et al.,, 2019).

Alongside particle aggregation, settling-driven gravitational instabilities
contribute to the early deposition of fine ash with similar outcomes (e.g.
grainsize bimodality, premature sedimentation of fine ash). These
instabilities generate downward-moving ash columns (fingers) which grow from
the base of the ash cloud (Figure 1) (Carazzo and Jellinek,, 2012; Manzella et
al.,, 2015; Scollo et al.,, 2017).

Figure 1: Gravitational instabilities observed at the base of a volcanic plume
during (A) the 2011 Gamalama eruption. (Credit: AP) and (B) the 2010 eruption
of Eyjafjallajökull (Manzella et al.,, 2015)

This phenomenon has the potential to enhance the sedimentation rate of fine
ash beyond the terminal fall velocity of individual particles, reducing the
residence time of fine ash in the atmosphere. Thus, a rigorous understanding
of these processes is important in order to build a comprehensive
parametrisation that can be included in dispersal models (Bonadonna et al.,,
2012; Durant,, 2015; Folch,, 2012; Scollo et al.,, 2010).

Settling-driven gravitational instabilities occur at the interface between an
upper, buoyant particle suspension, e.g., a volcanic ash cloud, and a lower,
denser fluid, e.g., the underlying atmosphere (Burns and Meiburg,, 2012;
Davarpanah Jazi and Wells,, 2016; Hoyal et al.,, 1999; Manzella et al.,,
2015). Whilst the initial density configuration is stable, particle settling
across the density interface creates a narrow unstable region called the
particle boundary layer (PBL) (Carazzo and Jellinek,, 2012). Once this attains
a critical thickness (Hoyal et al.,, 1999), a Rayleigh-Taylor-like instability
(Chandrasekhar,, 1961; Sharp,, 1984) can form on the interface between the PBL
and the lower layer, generating finger-like structures which propagate
downwards. A further critical condition for instability is that the particle
settling velocity $V_{s}$ must be smaller than the finger propagation velocity
$V_{f}$ (Carazzo and Jellinek,, 2012). Thus, the occurrence of the instability
enhances the sedimentation rate (Manzella et al.,, 2015; Scollo et al.,,
2017). Alternatively, if $V_{s}$ is greater than the propagation velocity of
fingers $V_{f}$, then particles settle individually before a PBL can form and
no instability occurs. Settling-driven gravitational instabilities have been
widely studied in laboratory experiments that simulate various natural
settings. Many experiments have considered an initial two-layer system, where
the particle suspension is initially separated from the underlying denser
layer by a removable horizontal barrier (Davarpanah Jazi and Wells,, 2016;
Fries et al.,, 2021; Harada et al.,, 2013; Hoyal et al.,, 1999; Manzella et
al.,, 2015; Scollo et al.,, 2017) whilst other experiments have involved
injection of the suspension into a density-stratified fluid at its neutral
buoyancy level (Carazzo and Jellinek,, 2012; Cardoso and Zarrebini,, 2001).
Similar instabilities can also be studied by allowing fine particles to
sediment through the free surface between water and air (Carey,, 1997;
Manville and Wilson,, 2004). Additionally, dimensional analysis has been used
to predict that the downward propagation velocity of the generated fingers is
given by (Carazzo and Jellinek,, 2012; Hoyal et al.,, 1999)

$V_{f}=\left[g\left(\frac{\rho_{PBL}-\rho_{f}}{\rho_{f}}\right)\right]^{\frac{2}{5}}\left[\frac{\pi
V_{s}\delta^{2}_{PBL}}{4}\right]^{\frac{1}{5}},$ (1)

where $\rho_{PBL}$ is the PBL bulk density, $\rho_{f}$ the underlying fluid
density, $g=9.81m.s^{-2}$ the gravitational acceleration and $\delta_{PBL}$
the PBL thickness, which by analogy with thermal convection (Turner,, 1973) is
taken to be (Hoyal et al.,, 1999)

$\delta_{PBL}=\left(\frac{Gr_{c}\nu^{2}}{g^{\prime}}\right)^{\frac{1}{3}},$
(2)

where $g^{\prime}=g\left(\rho_{PBL}-\rho_{f}\right)/\rho_{f}$, $\nu$ the
kinematic viscosity and $Gr_{c}$ a critical Grashof number (see Table S1 in
Supplementary Material for all acronyms and symbols used in this paper). The
reduced gravity $g^{\prime}$ describes the change in the gravitational
acceleration due to buoyancy forces. Continuing the analogy with thermal
convection, it has been proposed that $Gr_{c}=10^{3}$ (Hoyal et al.,, 1999),
although recent experimental observations suggests $Gr_{c}\approx 10^{4}$ may
be more accurate (Fries et al.,, 2021). Therefore, for known particle and
fluid properties, it is possible to predict whether collective settling will
occur and fingers subsequently form using the condition $V_{f}\sup V_{s}$
(Hoyal et al.,, 1999). According to this relation, the limit between
collective and individual settling occurs when $V_{f}=V_{s}$. However, the
transition is likely to be smooth, with a transitionary regime where both
fluid-like and particle-like settling occur at the same time, as suggested by
Harada et al., (2013).

For the initial two-layer configuration, Hoyal et al., (1999) also developed a
series of analytical mass-balance models predicting the average particle
concentration in the lower layer depending on whether the upper and lower
layers were convecting or not. In the case of a quiescent upper layer and a
convective lower layer (convection initiated by finger propagation), the
evolution of the mass of particles in the lower layer $M_{2}$ depends on the
balance between the mass flux of particles arriving from the upper layer
$\dot{M}_{in}$ and the mass flux of particle leaving by sedimentation
$\dot{M}_{out}$

$\frac{dM_{2}}{dt}=\dot{M}_{in}-\dot{M}_{out},$ (3)

where $t$ is time. Assuming that
$M_{2}\left(t\right)=Ah_{2}C_{2}\left(t\right)$, where $C_{2}\left(t\right)$
is the average particle concentration in the lower layer, Hoyal et al., (1999)
solved this equation using $\dot{M}_{in}=AV_{s}C_{1}\left(0\right)$,
$\dot{M}_{out}=AV_{s}C_{2}\left(t\right)$ and the initial condition
$C_{2}\left(0\right)=0$. Thus

$C_{2}\left(t\right)=C_{1}\left(0\right)\left[1-e^{-\frac{V_{s}}{h_{2}}t}\right],$
(4)

where $C_{1}\left(0\right)$ is the initial particle concentration in the upper
layer, $h_{2}$ the lower layer thickness and $A$ the horizontal cross section
of the tank.

Further studies of settling-driven gravitational instabilities have taken
theoretical approaches, such as using linear stability analyses to predict the
growth rate and characteristic wavelengths of the instability at very early
stages (Alsinan et al.,, 2017; Burns and Meiburg,, 2012; Yu et al.,, 2013).
Moreover, various numerical models simulating settling-driven gravitational
instability have also been developed (Burns and Meiburg,, 2014; Chou and
Shao,, 2016; Jacobs et al.,, 2013; Keck et al.,, 2021; Yamamoto et al.,,
2015). Most numerical approaches to this problem have used continuum-phase
models, where the coupling between particles and fluid is strong enough to
describe them as a single-phase (Burns and Meiburg,, 2014; Chou and Shao,,
2016; Yu et al.,, 2014). This Eulerian description is valid under the
assumptions of sufficiently small particles and a large enough number of
particles such that the drag and gravitational forces are in equilibrium. The
condition on the particle size can be quantified through the Stokes number
(Burgisser et al.,, 2005; Roche and Carazzo,, 2019)

$St=\frac{\rho_{p}D^{2}_{p}U}{18\mu L},$ (5)

where $\rho_{p}$ is the particle density, $D_{p}$ the particle diameter, $\mu$
the dynamic viscosity and $U$ and $L$ characteristic velocity and length
scales of the flow. For $St\inf 1$, the particles and fluid can be considered
coupled and, providing there are enough particles, the continuum approach is
valid.

The Eulerian description can be extended to multiple phases in order to
simulate their interaction (e.g., gas-liquid interaction) using adaptive mesh
refinements to resolve the phase interfaces (Jacobs et al.,, 2013). However,
for large particle diameters and small particle volume fractions, collective
behaviour no longer occurs and the continuum-phase method cannot be applied.
In this case, there is a need to explicitly model particle motion, taking the
drag force into consideration (Chou and Shao,, 2016; Yamamoto et al.,, 2015).
This paper presents an innovative method to implement a continuum model by
coupling the Lattice Boltzmann Method (LBM) with a low-diffusivity finite
difference (FD) scheme. This model takes advantage of the LBM capabilities to
simulate complex flows through uniform grids and thus, the ease of coupling
with finite difference methods. This hybrid model has been validated by
comparing the results with those from linear stability analysis and laboratory
experiments (Fries et al.,, 2021). The validated model then allows us to gain
new insights into the fundamental processes by exploring experimentally-
inaccessible regions of the parameter space. We first describe the general
framework and governing equations that describe settling-driven gravitational
instabilities, then the configuration of the validatory experiments to which
we apply the model. Next, we propose a numerical strategy involving a hybrid
model in order to solve the system of equations. We then go on to present the
linear stability analysis before finally describing and discussing the results
of our simulations.

## 2\. Methods

### 2.1. Problem formulation

The model consists of a three-way coupling between fluid momentum, fluid
density, and particle volume fraction, based on the assumption that the
particle suspension can be represented by a continuum concentration field.
Moreover, the particle drag force is in equilibrium with the gravitational
force such that the forcing term in the fluid momentum equation is equivalent
to a buoyant force term (Boussinesq approximation), which depends on the
particle volume fraction $\phi\left(\vec{x},t\right)$ (Burns and Meiburg,,
2014; Chou and Shao,, 2016; Yu et al.,, 2014). $\phi\left(\vec{x},t\right)$
satisfies the advection-diffusion-settling equation

$\frac{\partial\phi}{\partial
t}+\vec{\nabla}.[\phi(\vec{u}_{f}-V_{s}\vec{e}_{z})]=D_{c}\nabla^{2}\phi,$ (6)

where $\vec{u}_{f}\left(\vec{x},t\right)$ is the fluid velocity, $D_{c}$ the
particle diffusion coefficient, $\vec{e_{z}}$ the vertical unit vector and
$\vec{x}=\left(x,y,z\right)$ the position coordinate. The fluid is considered
incompressible meaning $\vec{\nabla}.\vec{u}_{f}=0$. Thus, equation 6 becomes

$\frac{\partial\phi}{\partial
t}+(\vec{u}_{f}-V_{s}\vec{e}_{z}).\vec{\nabla}\phi-\phi\vec{\nabla}.(V_{s}\vec{e}_{z})=D_{c}\nabla^{2}\phi.$
(7)

The particle settling velocity depends on the ambient fluid density $\rho$,
which in turn depends on any transported density-altering properties, such as
temperature or the concentration of a chemical species, e.g., the sugar in our
validatory experiments (Fries et al.,, 2021). We incorporate the effect of a
single density-altering property on the fluid density through a classical
advection-diffusion equation

$\frac{\partial\rho\left(\rho_{0},S\right)}{\partial
t}+\vec{u}_{f}.\vec{\nabla}\rho\left(\rho_{0},S\right)=D_{s}\nabla^{2}\rho\left(\rho_{0},S\right),$
(8)

where $\rho_{0}$ is a reference density of the carrier fluid, $S$ the density-
altering quantity (temperature or concentration), and $D$ the associated
diffusion coefficient. Additionally, under the Boussinesq approximation, we
assume that the density depends linearly on $S$. The fluid momentum is
modelled with the incompressible Navier-Stokes momentum equation

$\frac{\partial\vec{u}_{f}}{\partial
t}+(\vec{u}_{f}.\vec{\nabla})\vec{u}_{f}=-\frac{1}{\rho_{0}}\vec{\nabla}p+\nu\nabla^{2}\vec{u}_{f}+\vec{F},$
(9)

where $p$ is the pressure and $\vec{F}$ the buoyant body force term. We
complete the system of equations by taking this force term to be a function of
$\phi$ and $\rho$

$\vec{F}=\left[\left(\frac{\rho_{p}-\rho_{0}}{\rho_{0}}\right)\phi+\left(\frac{\rho}{\rho_{0}}-1\right)(1-\phi)\right]\vec{g}.$
(10)

The system of equations presented so far assumes that all particles are of
uniform size. In order to generalise to systems with polydisperse particle
size distributions, we consider $N$ different particle concentration fields
$\phi_{i}$, where each one is associated with a different size class and
individually satisfies equation 7. Furthermore, the body force term becomes

$\vec{F}=\left[\left(\frac{\rho_{p}-\rho_{0}}{\rho_{0}}\right)\phi_{tot}+\left(\frac{\rho}{\rho_{0}}-1\right)(1-\phi_{tot})\right]\vec{g}$
(11)

where

$\phi_{tot}=\sum_{i=1}^{N}\phi_{i}.$ (12)

### 2.2. Flow configuration and experiment description

Full details of the validatory laboratory experiments can be found in Fries et
al., (2021) but we summarise the essential details here. The experiments are
performed in a configuration identical to that of Manzella et al., (2015) and
Scollo et al., (2017) (Figure 2) and consist of a water tank divided into two
horizontal layers, initially separated by a removable barrier. The upper layer
is an initially mixed particle suspension, which represents the ash cloud, and
the lower layer is a dense sugar solution, analogue to the underlying
atmosphere. The particles are spherical glass beads with a median diameter of
$41.5\pm 0.5\,\mu m$ (measured using laser diffraction with a Bettersizer S3
Plus) and a density $\rho_{p}$ of $2519.4\pm 0.09\,kg/m^{3}$ (measured using
helium pycnometry UltraPyc 1200e), and are sufficiently small to be well-
coupled with the fluid, whilst the initial particle concentration
$C_{1}\left(0\right)$ of the upper layer is varied from $1$ to $10\,g/l$ (see
Table 1 for the conversion to particle volume fraction $\phi_{0}$). The lower
layer density is kept constant at $\rho_{f}=1008.4\,kg/m^{3}$ (corresponding
to a sugar concentration of $S_{0}=35\,g/l$), always ensuring an initially
stable density configuration.

At the beginning of an experiment, the barrier separating the two layers is
removed, allowing particle settling through the interface. A PBL subsequently
forms and finger formation is initiated. Experiments are illuminated from the
side of the water tank with a planar laser and recorded with a high-contrast
camera. We measure the vertical finger velocity by tracking the progression of
the finger front with time. Additionally, Planar Laser Induced Fluorescence
(PLIF) (Crimaldi,, 2008; Koochesfahani,, 1984) and particle imaging are used
to quantify the spatial distribution of the fluid phase density and particle
concentration.

Figure 2: Experimental setup used by Fries et al., (2021) and the initial density profiles associated with the contributions from particles (blue dashed) and sugar (red dotted), as well as the bulk density (black solid). The density of fresh water is given by $\rho_{0}$ Initial particle concentration $(g/l)$ |  Volume fraction |  Particle diameter $(\mu m)$ | $z_{H}\,(mm)$ | $z_{\phi}\,(mm)$ | $z_{s}\,(mm)$
---|---|---|---|---|---
1 |  $3.97\times 10^{-4}$ |  40 | 12.11 | 2.59 | 0.67
2 |  $7.94\times 10^{-4}$ |  40 | 7.37 | 2.20 | 0.63
3 |  $1.19\times 10^{-3}$ |  40 | 5.47 | 1.99 | 0.61
4 |  $1.59\times 10^{-3}$ |  40 | 5.01 | 1.99 | 0.61
5 |  $1.98\times 10^{-3}$ |  40 | 4.20 | 1.84 | 0.60
6 |  $2.38\times 10^{-3}$ |  40 | 3.45 | 1.72 | 0.59
7 |  $2.78\times 10^{-3}$ |  40 | 3.25 | 1.70 | 0.60
8 |  $3.18\times 10^{-3}$ |  40 | 3.17 | 1.72 | 0.61
9 |  $3.57\times 10^{-3}$ |  40 | 3.01 | 1.71 | 0.62
10 |  $3.97\times 10^{-3}$ |  40 | 2.54 | 1.58 | 0.61
3 |  $1.19\times 10^{-3}$ |  25 | - | - | -
3 |  $1.19\times 10^{-3}$ |  55 | - | - | -
3 |  $1.19\times 10^{-3}$ |  70 | - | - | -
3 |  $1.19\times 10^{-3}$ |  85 | - | - | -
3 |  $1.19\times 10^{-3}$ |  100 | - | - | -
3 |  $1.19\times 10^{-3}$ |  115 | - | - | -
3 |  $1.19\times 10^{-3}$ |  130 | - | - | -
9 |  $3.57\times 10^{-3}$ |  25 | - | - | -
9 |  $3.57\times 10^{-3}$ |  55 | - | - | -
9 |  $3.57\times 10^{-3}$ |  70 | - | - | -
9 |  $3.57\times 10^{-3}$ |  85 | - | - | -
9 |  $3.57\times 10^{-3}$ |  100 | - | - | -
9 |  $3.57\times 10^{-3}$ |  115 | - | - | -
9 |  $3.57\times 10^{-3}$ |  130 | - | - | -
9 |  $3.57\times 10^{-3}$ |  145 | - | - | -
9 |  $3.57\times 10^{-3}$ |  160 | - | - | -
Table 1: List of simulations performed. All the simulations have been
performed using an initial lower layer fluid density of $1008.4\,kg/m^{3}$.
$z_{H}$, $z_{\phi}$ and $z_{s}$ are parameters used in the linear stability
analysis (LSA) in order to describe the different base states associated with
the particle and sugar profiles in equations 39 and 40. The LSA has been
performed only for a constant particle size of $40\,\mu m$ in order to study
the effect of the particle volume fraction.

####  2.2.1 Application to flow configuration

We apply the general system of equations presented in section 2.1 to the
configuration of the validatory experiments. The particles are spherical and
sufficiently small that their terminal settling velocity in water is given by
the Stokes velocity (Stokes,, 1851)

$V_{s}=\frac{D^{2}_{p}g\left[\rho_{p}-\rho\left(S\right)\right]}{18\mu},$ (13)

where $S$ is the sugar concentration and $\rho=\rho_{0}\left(1+\alpha
S\right)$, with $\alpha$ the sugar expansion coefficient.

We simulate the solid walls of the tank around our domain with a no-slip
boundary condition for the fluid velocity. Neumann boundary conditions are
employed for $\phi$ and $\rho$ to avoid any flux of particles or sugar across
the walls. Thus we impose

$\frac{\partial\phi}{\partial x}=0,$ (14)

and

$\frac{\partial\rho}{\partial x}=0,$ (15)

on the wall nodes. Furthermore, we define the following initial states for
$\phi$ and $S$

$\phi\left(\vec{x},t=0\right)=\left\\{\begin{array}[]{rcl}0&\mbox{for}&z<H_{0}\\\
\phi_{0}&\mbox{for}&z>H_{0}\end{array}\right.$ (16)

and

$S\left(\vec{x},t=0\right)=\left\\{\begin{array}[]{rcl}S_{0}&\mbox{for}&z<H_{0}\\\
0&\mbox{for}&z>H_{0}\end{array}\right.$ (17)

where $\phi_{0}$ and $S_{0}$ are the initial particle volume fraction in the
upper layer and initial sugar concentration in the lower layer, respectively,
and $H_{0}=0.25\,m$ the initial height of the interface ($z=0$ corresponds to
the base of the tank). We also add a small perturbation to the particle volume
fraction field in order to initiate the instability. Finally, the system is
initially stationary so $\vec{u_{f}}\left(\vec{x},t=0\right)=0$.

### 2.3. Numerical methods

The model is implemented using a hybrid strategy where a LBM solves the fluid
motion and is coupled with finite difference schemes that solve the advection-
diffusion equations for $\phi$ and $S$.

####  2.3.1 Fluid motion

The LBM is an efficient alternative to conventional Computational Fluid
Dynamics (CFD) methods that explicitly solve the Navier-Stokes equations at
each node of a discretised domain (Chen and Doolen,, 2010; He and Luo,, 1997).
It is a well-established approach for simulating complex flows, including
multiphase fluids (Leclaire et al.,, 2017) and thermal and buoyancy effects
(Noriega et al.,, 2013; Parmigiani et al.,, 2009). The LBM originates from the
kinetic theory of gases and provides a description of gas dynamics at the
mesoscopic scale. This scale exists between the microscopic, which describes
molecular dynamics, and the macroscopic, which gives a continuum description
of the system with variables such as density and velocity. Thereby, the
mesoscopic scale considers a probability distribution function of molecules
described by the Lattice Boltzmann equation. This model reduces the process to
two main steps: streaming (i.e., displacement of populations between
consecutive calculation nodes), and collision (i.e., interaction of
populations on a node). The Bhatnagar-Gross-Krook (BGK) model (Bhatnagar et
al.,, 1954) provides a simple collision process based on a fundamental
property given by kinetic theory which describes gas motion as a perturbation
around the equilibrium state. Then, the LBM-BGK model solves, for the particle
population $f_{i}$, which are a discrete representation of the probability
distribution function, the equation

$f_{i}\left(\vec{x}+\vec{c_{i}}\delta t,t+\delta t\right)=-\frac{\delta
t}{\tau}\left(f_{i}\left(\vec{x},t\right)-f^{eq}_{i}\right),$ (18)

where $\delta t$ is the time step, $f^{eq}_{i}\left(\rho,\vec{u_{f}}\right)$
the equilibrium distribution function, $\tau$ the relaxation time associated
with the flow viscosity and $\vec{c_{i}}$ the local particle velocity. The LBM
is applied to specific types of lattices that describe how the populations
move through the calculation nodes (Kruger et al.,, 2017). These types of
lattice are commonly summarized in the form $DrQm$ where $r$ denotes the
dimension of the system and $m$ the number of directions in which populations
can propagate.

Figure 3 shows the scheme D3Q19 used for our 3D simulations and the associated
set of local velocities.

Figure 3: Depiction of the D3Q19 lattice. The red arrows show the different
possible directions of propagation. The associated local velocities are
summarised in the velocity set $\vec{c_{i}}$.

The macroscopic fluid state is described through the usual macroscopic
variables such as density, velocity and kinematic viscosity. These variables
are related to the moments of the populations $f_{i}$ through

$\rho=\sum_{i}f_{i},$ (19)

and

$\rho\vec{u_{f}}=\sum_{i}f_{i}\vec{c_{i}},$ (20)

whilst the kinematic viscosity controls the relaxation to equilibrium through
the relaxation time

$\tau=\frac{\nu}{c^{2}_{s}}+\frac{\delta t}{2},$ (21)

The variable $c_{s}^{2}$ is commonly called the speed of sound and is equal to
$(1/3)(\delta x^{2}/\delta t^{2})$ where $\delta x$ is the spatial step.
However, the classical LBM-BGK model described above does not take into
account any forcing term. One way to include forcing is to rewrite equation 18
as

$f_{i}\left(\vec{x}+\vec{c_{i}}\delta t,t+\delta t\right)=-\frac{\delta
t}{\tau}\left(f_{i}\left(\vec{x},t\right)-f^{eq}_{i}\right)+F_{i}\delta t,$
(22)

where $F_{i}$ is the forcing term, which can be expressed by a power series in
the local particle velocities, and the equilibrium distribution is now given
by $f^{eq}_{i}=f^{eq}_{i}\left(\rho,\vec{u^{*}_{f}}\right)$, where
$\rho\vec{u^{*}_{f}}=\sum_{i}f_{i}\vec{c_{i}}+b\vec{F}\delta t$. The
determination of the coefficient $b$, as well as the power series expansion of
$F_{i}$ are described by Guo et al., (2002). Finally, no-slip boundary
conditions in the LBM, to simulate walls for example, can be implemented using
the classical $bounce-back$ boundary condition (Kruger et al.,, 2017) where
the populations arriving on a wall node during the streaming step are simply
reflected back to their previous nodes.

####  2.3.2 Transport of particles and other density-altering quantities

The particles and other density-altering quantities are described by continuum
fields that follow an advection-diffusion law coupled with the fluid motion as
simulated with the LBM. The numerical solution of the advection equation is
particularly challenging for methods which, like ours, are Eulerian (i.e.,
mesh-based). Indeed, such methods exhibit numerical diffusion which may
strongly reduce model accuracy and, in some cases, even exceed the amplitude
of the actual, physical diffusion term. The lack of physical diffusion in our
problem and the presence of sharp interfaces restrict our ability to solve the
advection equations with the LBM. In fact, the advection-diffusion equation
can be solved by the LBM with a BGK approach in analogous fashion to the fluid
motion by modifying the equilibrium distribution and the relaxation time to
depend on the diffusion coefficient $D$ rather than $\nu$

$\tau=\frac{D}{c^{2}_{s}}+\frac{\delta t}{2},$ (23)

However, a stability condition for a LBM-BGK algorithm is $\tau/\delta t=1/2$.
Thus, since the problem is convection dominated, the low diffusion coefficient
($D\ll 1$) drives the model towards the stability limit, introducing strong
numerical errors near sharp concentration gradients (Hosseini et al.,, 2017).
For this reason, we solve the advection term using two finite-difference
schemes which are selected depending on the required accuracy: the classical
first-order upwind finite difference and the third-order Weighted Essentially
Non Oscillatory (WENO) finite difference scheme (Jiang and Shu,, 1996; Liu et
al.,, 1994).

Coupling the LBM with an upwind finite difference scheme allows us to avoid
the stability problem. First-order FD schemes however, still suffer from the
problem of numerical diffusion due to the truncation error associated with
terminating the Taylor expansion after the first spatial derivative. The
induced numerical error NE for the convective term in the advection-diffusion
equation is given by

$NE\sim u\frac{\delta x}{2}\frac{\partial^{2}\phi}{\partial x^{2}},$ (24)

where $u$ is the transport velocity. $NE$ acts like an additional diffusion
term because of the presence of the second-order derivative. The numerical
diffusion associated with the solution of $S$ is negligible due to the low
fluid velocity and consequently the use of the first order FD scheme does not
significantly affect the accuracy. However, in the solution of $\phi$, which
includes an additional velocity contribution due to the settling, the
truncation error associated with the first-order scheme becomes non-
negligible. Whilst decreasing $\delta x$ would reduce numerical diffusion, we
would require an unpractically small value in order to get a sufficiently
accurate solution. Additionally, simply increasing the order of the scheme
introduces dispersion (spurious oscillations) near regions of high gradient,
according to the Godunov theorem (Godunov,, 1954, 1959). Therefore, we choose
here to implement the low diffusive WENO procedure for the solution of $\phi$,
thus achieving a stable and high-resolution scheme without dispersion.

Further information on how we discretise the convective term in the advection-
diffusion equation using the first order upwind and the third order WENO
finite difference schemes is detailed in Section 1 of the Supplementary
material.

####  2.3.3 Numerical implementation

Our model is implemented using $Palabos$ (Parallel Lattice Boltzmann Solver),
a Computational Fluid Dynamics (CFD) solver based on the Lattice Boltzmann
Method and developed by the Scientific Parallel Computing group of the
Computer Science Department, University of Geneva (Latt et al.,, 2020).
$Palabos$ is designed to perform calculations on massively parallel computers,
thus allowing very small spatial resolutions in order to accurately simulate
the finger dynamics.

## 3\. Linear stability analysis

In order to validate our model, we compare the early-time simulated behaviour
against predictions from linear stability analysis (LSA). LSA is applied to
the onset of the physical instability at the interface between layers of
different particle concentration. It involves defining a field equation-
satisfying base state for each of the unknown fields in a problem and then
applying an infinitesimally small perturbation to each of these fields. The
equations are then expanded to linear order in the perturbation, with higher
order terms assumed to be negligible. By assigning the perturbation to have
the form of a complex waveform, the system of equations reduces to an
eigenvalue problem, which can be solved to determine which wavelengths will
grow or decay (Chandrasekhar,, 1961). In this section, we assume that the
system is invariant under translation in the x-y plane, thus reducing the
analysis to a 2D problem. We strongly follow the procedure described by (Burns
and Meiburg,, 2012) in order to solve our problem.

### 3.1. Nondimensionalisation

We nondimensionalise our system of equations by defining

$l^{c}=\left(\frac{\nu^{2}}{g}\right)^{\frac{1}{3}},$ (25)

$t^{c}=\left(\frac{\nu}{g^{2}}\right)^{\frac{1}{3}},$ (26)

and

$p^{c}=\rho_{0}\left(\nu g\right)^{\frac{2}{3}},$ (27)

where $l^{c}$, $t^{c}$ and $p^{c}$ are characteristic quantities. We also
define the dimensionless parameters

$S^{*}=\alpha S_{0},$ (28)

$\phi^{*}=\phi_{0},$ (29)

$Fr=\frac{1}{t^{c}}\sqrt{\frac{l^{c}}{g}},$ (30)

and

$Sc_{i}=\frac{\nu}{D_{i}},$ (31)

noting that $Fr$ is a Froude number and $Sc_{i}$ are Schmidt numbers.
Furthermore, the stream function $\psi$ is defined such that
$\vec{u_{f}}=\left(\partial\psi/\partial z,-\partial\psi/\partial x\right)$
and the vorticity as $\vec{\omega}=\vec{\nabla}\times\vec{u_{f}}$. Then,
applying the characteristic quantities to the vorticity formulation and
equations (7-9), we obtain the dimensionless system (for the rest of the
analysis, all the symbols used represent dimensionless quantities)

$\omega=-\nabla^{2}\psi,$ (32)

$\frac{\partial\omega}{\partial
t}+\left(\vec{u_{f}}.\vec{\nabla}\right)\omega=\nabla^{2}\omega+\frac{\partial\phi}{\partial
y}\frac{\phi^{*}}{Fr^{2}}\left[SS^{*}-\left(\frac{\rho_{p}-\rho_{0}}{\rho_{0}}\right)\right]-\frac{\partial
S}{\partial y}\frac{S^{*}}{Fr^{2}}\left(1-\phi\phi^{*}\right),$ (33)

$\frac{\partial S}{\partial
t}+\vec{u_{f}}.\vec{\nabla}S=\frac{1}{Sc_{s}}\nabla^{2}S,$ (34)

and

$\frac{\partial\phi}{\partial
t}+\left(\vec{u_{f}}-V_{s}\vec{e_{z}}\right).\vec{\nabla}\phi=\frac{1}{Sc_{c}}\nabla^{2}\phi.$
(35)

Note that here we have neglected the term
$-\phi\vec{\nabla}.\left(V_{s}\vec{e_{z}}\right)$ in equation 7 assuming the
fluid density variation across the interface is sufficiently small that it
does not affect the particle settling velocity.

### 3.2. Variable expansion and eigenvalue problem

We linearise the system of equations by expanding each variable in terms of a
base state and a perturbation

$\varphi\left(y,z,t\right)=\bar{\varphi}\left(z\right)+\varphi^{\prime}\left(y,z,t\right),$
(36)

where $\varphi\left(y,z,t\right)=\\{\psi,\omega,\phi,S\\}$,
$\bar{\varphi}\left(z\right)=\\{\bar{\psi},\bar{\omega},\bar{\phi},\bar{S}\\}$
the associated base state and
$\varphi^{\prime}\left(y,z,t\right)=\\{\psi^{\prime},\omega^{\prime},\phi^{\prime},S^{\prime}\\}$
the perturbation. We choose the following base states

$\bar{\psi}=0,$ (37)

$\bar{\omega}=0,$ (38)

$\bar{\phi}\left(z,t\right)=\frac{1}{2}\left[1+erf\left(\frac{z}{z_{\phi}}\right)\right],$
(39)

$\bar{S}\left(z,t\right)=\frac{1}{2}\left[1-erf\left(\frac{z-V_{s}t}{z_{S}}\right)\right],$
(40)

where $z_{\phi}\left(T\right)$ and $z_{S}\left(T\right)$ are coefficients
fitted in order to have similar base states to the profiles observed in the
simulations prior to the onset of the instability which starts growing at the
time $T$.

Solutions for the perturbation are assumed to have the form of normal modes

$\varphi^{\prime}\left(y,z,t\right)=\hat{\varphi}\left(z\right)exp\left(iky+\sigma
t\right),$ (41)

where $\hat{\varphi}\left(z\right)$ is the perturbation amplitude, $k$ the
wavenumber and $\sigma$ the instability growth rate. The linearised system of
equations is then formulated in matrix form so that the problem is reduced to
the eigenvalue problem $K\vec{x}=\sigma W\vec{x}$ where

$\vec{x}=\left(\begin{array}[]{c}\hat{\psi}\left(z\right)\\\
\hat{\omega}\left(z\right)\\\ \hat{S}\left(z\right)\\\
\hat{\phi}\left(z\right)\end{array}\right),$ (42)

and, in a reference frame moving downward at $V_{s}$, the matrices $K$ and $W$
are given by

$K=\left(\begin{array}[]{c c c c}M&I&0&0\\\
0&M-V_{s}D_{z}&-ik\frac{S^{*}}{Fr^{2}}\left(1-\bar{\phi}\phi^{*}\right)I&ik\frac{\phi^{*}}{Fr^{2}}\left[SS^{*}-\left(\frac{\rho_{p}-\rho_{0}}{\rho_{0}}\right)\right]I\\\
ik\frac{d\bar{S}}{dz}I&0&\frac{1}{Sc_{S}}M-V_{s}D_{z}&0\\\
ik\frac{d\bar{\phi}}{dz}I&0&0&\frac{1}{Sc_{c}}M\end{array}\right),$ (43)

and

$W=\left(\begin{array}[]{c c c c}0&0&0&0\\\ 0&I&0&0\\\ 0&0&I&0\\\
0&0&0&I\end{array}\right),$ (44)

where $D_{z}=\partial/\partial z$, $M=-k^{2}+D^{2}_{z}$ and $I$ is the
identity operator.

The eigenvalues $\sigma$ determine the stability of the system:

  * •

If all the eigenvalues have negative real parts, the system remains stable

  * •

If at least one eigenvalue has a positive real part, the system is unstable.

In order to solve the eigenvalue problem, the spatial derivatives are
discretised using the linear rational collocation method with a grid
transformation allowing a fine resolution around narrow interfacial regions
(Baltensperger and Berrut,, 2001; Berrut and Mittelmann,, 2004).

The key result of the LSA is the dispersion relation between $\sigma$ and $k$.
Figure 4 presents the growth rate as a function of the wavenumber, for
different initial particle volume fractions. The parameters for the different
base states used to produce these curves are summarised in the Table 1. We use
this result in Section 4.1 in order to compare the predictions of the LSA with
the results of our numerical model.

Figure 4: Dispersion relation obtained from LSA for several initial particle
concentrations.

## 4\. Results

We validate our numerical model by comparing the results with predictions from
LSA and experimental observations. The LSA predicts the growth rates of
different perturbation wavenumbers during the very early stage of the
instability, which can be compared with the spectrum of wavenumbers present in
the particle concentration interface in the numerical model. Additionally, the
experiments of Fries et al., (2021) employ imaging techniques to measure
quantities, such as the particle concentration field and finger velocity, at
times beyond the linear regime. Finally, our results are compared with some
results of previous analyses on settling-driven gravitational instabilities
(Carazzo and Jellinek,, 2012; Hoyal et al.,, 1999).

### 4.1. Comparison of model results with predictions from linear stability
analysis

In order to compare our 3D simulations with the 2D linear stability analysis,
we consider just the central plane of the simulation domain, i.e., a slice in
the $\left(y,z\right)$ plane located at $x=l_{x}/2$ ($l_{x}$ being the tank
depth) (Figure 2). We define the front of the particle field to be the lowest
position where $\phi=\phi_{0}/2$ and also define $H\left(y\right)$ to be the
separation between $z=0$ and this front. Our study has shown that the front
position is only weakly affected when using other possible thresholds, i.e.,
$\phi_{0}/10$ or $\phi_{0}/5$ (relative change $\leavevmode\nobreak\ 3\%$).
Figure 5A shows an example of a space-time diagram showing the evolution of
$H\left(y\right)$ through time. Furthermore, by calculating the Fourier
transform $\tilde{H}\left(k,t\right)$ of $H\left(y\right)$ at different times,
we can identify different dominant wavenumbers and their associated amplitudes
as shown in the space diagram of the power spectral density (PSD)
$\Gamma_{H}\left(k,t\right)=\left(1/\left(k_{S}L_{S}\right)\right)|\tilde{H}\left(k,t\right)|^{2}$
(Figure 5B), where $k_{S}$ is the sampling wavenumber and $L_{S}$ the number
of samples. We extract the dominant mode and its associated growth rate from
$\Gamma_{H}\left(k,t\right)$ and compare the results with the predicted growth
rates from LSA. We apply this analysis during a period when the amplitude
$|\tilde{H}\left(k,t\right)|$ of any given mode does not exceed $40\%$ of its
wavelength, thus ensuring we are still in the linear regime (Lewis,, 1950).

During the linear regime, we can assume that the growth of the spectral
amplitude can be described as (Völtz et al.,, 2001)

$|\tilde{H}\left(k,t\right)|=|\tilde{H}_{i}\left(k\right)|exp\left(\sigma_{sim}\left(k\right)t\right),$
(45)

with $|\tilde{H}_{i}\left(k\right)|$ the initial amplitude and
$\sigma_{sim}\left(k\right)$ the instability growth rate as determined from
the simulations. Thus, the PSD can be expressed as

$\Gamma_{H}\left(k,t\right)=\Gamma_{H_{i}}exp\left(2\sigma_{sim}\left(k\right)t\right),$
(46)

where $\Gamma_{H_{i}}$ is the initial spectral density.

Figure 5: (A) Space-time diagram of the particle front height
$H\left(y,t\right)$. (B) Evolution of the power spectral density of the
particle interface over time. Initial particle volume fraction:
$\phi_{0}=3.97\times 10^{-4}$.

At each time step, we extract the PSD and the wavenumber $k_{sim}$ associated
with the dominant mode as shown in Figure 6. However, we observe that the
dominant mode remains at the same wavenumber during instability growth except
for three cases ($\phi_{0}=1.59\times 10^{-3}$, $2.38\times 10^{-3}$ and
$3.97\times 10^{-3}$) where we observed that the dominant mode changed its
position in the spectral space. For these simulations only, we have a set of
several wavenumbers $k_{sim,i}$, ($i=1,2,3$) associated with the dominant
mode. With the computed PSD of the dominant mode as a function of time
$\Gamma_{H}\left(k_{sim},t\right)$, we apply our exponential fitting (equation
46) to determine the growth rate $\sigma_{sim,i}$ (Figure 6B). For the
simulations which resulted in several values of $k_{sim,i}$ for the dominant
mode, we measured the growth rates of each mode $\sigma_{sim,i}$ and found
identical values, up to a precision of $5\%$. Additionally, for each
simulation, we find the time $T$ when the instability starts growing, easily
identified as the time at which the modal wavenumber becomes non-zero (e.g.,
in Figure 6A this is at approximately 6 s). At this time, we extract the
associated vertical profiles of particle and sugar concentration which are
used to find the coefficients $z_{\phi}(T)$ and $z_{S}(T)$ (equations 39-40)
and thus determine the base states of $\bar{\phi}(z,T)$ and $\bar{S}(z,T)$
(Figure 7). We then perform the LSA for each $\phi_{0}$, using the appropriate
base states, and obtain a dispersion relation $\sigma=f(k)$. Using this
relation, we predict the different growth rates $\sigma=\sigma_{LSA,i}$
associated with $k=k_{sim,i}$ and we compare with $\sigma(k_{sim,i})$ as
measured in our simulations. Figure 8 shows the comparison between
$\sigma_{sim}$ (black dots) and $\sigma_{LSA,i}$ (red triangles), as predicted
from the LSA, for the dominant wave mode. For the cases including a moving
dominant mode, we plotted the growth rates associated with the different
measured wavenumbers. We see that the dependence of the largest value of
$\sigma_{LSA,i}$ on the initial particle concentration is in good agreement
with the simulated growth rate.

Figure 6: (A) Example of dominant wavenumbers extracted from the maximum of
the PSD. Initial particle volume fraction $\phi_{0}=7.94\times 10^{-4}$. (B)
Exponential fitting to the temporal evolution of the PSD for the first maximum
in (A), $k_{sim,1}=0.517\,mm^{-1}$. Figure 7: Example of LSA base states
extracted from the simulations for $\phi_{0}=3.97\times 10^{-4}$. Dots:
profiles extracted from the simulations at $T=9.55s$ (start of the instability
growth). Dotted lines: Fit with equations 39 and 40. i.e., base states used
for the LSA. Blue: particle volume fraction. Red: Sugar concentration. Figure
8: Comparison of the instability growth rate measured in the simulations
(black circles) and that predicted by the linear stability analysis (red
triangles).

### 4.2. Comparison with experimental investigations

Figures 9A and 9B show a qualitative comparison between snapshots taken from
experiments (Fries et al.,, 2021) and simulations (slice in the numerical
domain). First, we note that our model is able to qualitatively reproduce the
shape and size of fingers, especially their fronts where we observe the
formation of lobes and eddies due to the Kelvin-Helmholtz instability (Chou
and Shao,, 2016). Second, we provide a quantitative validation of the non-
linear regime by comparing our model with experiments, through measurements of
the PBL thickness and the vertical finger velocity as functions of the
particle volume fraction and size.

Figure 9: Settling-driven gravitational instabilities observed 19.5 s after
the barrier removal (A) in the laboratory (Fries et al.,, 2021) and (B) in
numerical simulations. Particle size: $40\,\mu m$ and initial volume fraction
: $\phi_{0}=2.78\times 10^{-3}$.

####  4.2.1 Characterisation of the PBL and effect of the initial particle
volume fraction on the finger velocity

The bulk density profile $\rho_{blk}$, derived from the contributions of the
particle concentration and sugar profiles, is given everywhere by the relation

$\rho_{blk}=\phi\rho_{p}+\left(1-\phi\right)\rho_{f}.$ (47)

Figure 10 shows the profiles of $\phi$, $\rho_{f}$ and $\rho_{blk}$ in the
numerical simulations as well as in the experiments 8 seconds after the
barrier removal for the same initial conditions ($\phi_{0}=3.18\times 10^{-3}$
). It can be seen that, in both the model and the experiments, there is an
increase of the bulk density below the initial interface, owing to the
particle front moving downwards. This zone of excess density corresponds to
the unstable PBL from which instabilities occur, generating fingers. To
calculate the finger velocity using the same method as in experiments, we
extract slices from the 3D numerical domain and manually track the fronts of
several fingers (6 to 15 fingers) from when they become fully developed until
just before they become too diluted. For each simulation with different volume
fraction, we then average the velocity of all tracked fingers and the
uncertainty is the standard deviation associated with each set of fingers used
for the measurements. Figure 11A shows the average finger velocity $V_{f}$ as
a function of $\phi_{0}$, for both experiments (Fries et al.,, 2021) and
simulations. Our simulation results are in good agreement with the
experimental measurements and highlight that the increase of $V_{f}$ with
$\phi_{0}$ is non-linear.

Figure 10: Density profile after 8 seconds for experiments (left) (Fries et
al.,, 2021) and simulations (right) with $\phi_{0}=3.18\times 10^{-3}$ and a
particle size of $40\,\mu m$. Figure 11: (A) Average finger speed ($V_{f}$) as
a function of the initial volume fraction ($\phi_{0}$) for a particle diameter
of $40\,\mu m$. Red and black dotted lines show the best fits to the
experiments (Fries et al.,, 2021) and simulations, respectively, using
equation 1 with $Gr_{c}$ as the fit parameter. For the simulations, we find
$Gr_{c}=1.2\pm 0.4\times 10^{4}$ whilst for the experiments $Gr_{c}=1.9\pm
0.7\times 10^{4}$. (B) Average finger speed ($V_{f}$) as a function of the
initial particle diameter ($D_{p}$), for two different particle volume
fractions. The green line is the Stokes velocity for individual particles. The
black dotted lines show the best fits to the simulations using equation 1 with
$Gr_{c}$ as the fit parameter. For $\phi_{0}=1.19\times 10^{-3}$, the best fit
gives $Gr_{c}=7.6\times 10^{3}$ and no fingers are observed to form for
particle sizes higher than $115\,\mu m$. For $\phi_{0}=3.57\times 10^{-3}$,
the best fit gives $Gr_{c}=2.7\times 10^{4}$ and no fingers are observed to
form for particle sizes higher than $145\,\mu m$. In the two plots, the blue
dashed line shows equation 1 using $Gr_{c}=10^{3}$ from the analogy with
thermal convection (Hoyal et al.,, 1999).

By analogy with thermal convection, it has previously been assumed that
$Gr_{c}=10^{3}$ (Hoyal et al.,, 1999), but this is only an order of magnitude
estimate and its application to settling-driven gravitational instabilities
remains uncertain (Fries et al.,, 2021). Figure 11A shows good agreement
between the simulations and equation 1 for a fitted $Gr_{c}$ of $1.2\pm
0.4\times 10^{4}$ ($R^{2}=0.92$), which is an order of magnitude higher than
the value previously assumed (Carazzo and Jellinek,, 2012; Hoyal et al.,,
1999). This agrees reasonably with the experiments, where the best fit is
obtained for $Gr_{c}=1.9\pm 0.7\times 10^{4}$ ($R^{2}=0.75$), but the
experimental results show more scatter. However, neither of these fits have
completely satisfactory values of $R^{2}$. We therefore further investigate
the applicability of equation 1 by examining the dependence of $V_{f}$ on
$\phi_{0}$, assuming a power law of the form $V_{f}\propto\phi_{0}^{q}$.
According to equation 1, $q=4/15\approx 0.27$. However, from the experiments,
we obtain $q=0.50\pm 0.16$ (with $R^{2}=0.95$) while for our simulations
$q=0.37\pm 0.08$ (with $R^{2}=0.98$). Here $Gr_{c}$, $q$ and their associated
uncertainties have been calculated accounting for the uncertainty on $V_{f}$
with the SciPy (Python-based ecosystem) procedure
$scipy.optimize.curve_{f}it$.

####  4.2.2 Effect of particle size on the finger vertical velocity

Since gravitational instabilities cause particles to sediment faster than
their settling velocity, it is of interest to explore the transition from
collective to individual settling, since this has implications for which grain
sizes may prematurely sediment from a volcanic cloud (Scollo et al.,, 2017).
Figure 11B shows the effect of particle size on the finger velocity as
measured from the model, for two different initial volume fractions, in the
experiments configuration (i.e. in the tank filled with water). We clearly
observe two regimes:

  * •

For particle sizes less than or equal to $115\,\mu m$ (for
$\phi_{0}=1.19\times 10^{-3}$) and $145\,\mu m$ (for $\phi_{0}=3.57\times
10^{-3}$), we observe fingers, with the finger velocity increasing with
particle size.

  * •

For greater particle sizes, no fingers are observed to form.

From our simulations, we constrain the transition between the two regimes to
occur at a critical particle diameters around $115\,\mu m$ and $145\,\mu m$
respectively for $\phi_{0}=1.19\times 10^{-3}$ and $\phi_{0}=3.57\times
10^{-3}$. We also note that this size range corresponds to the particle size
at which the Stokes velocity exceeds the predicted finger velocity. This
result agrees with the experimental observations of Scollo et al., (2017), who
observed that no fingers form for particles with diameter larger than $\sim
125\,\mu m$ with in initial particle volume fraction of $\phi_{0}=1.19^{-3}$.
We also compare the dependence of $V_{f}$ on the particle diameter with that
predicted by equation 1 and find a best fit for $Gr_{c}=7.6\pm 3.6\times
10^{3}$ (with $R^{2}=0.91$) and $Gr_{c}=2.7\pm 0.8\times 10^{4}$ (with
$R^{2}=0.87$) respectively for the two initial volume fractions (Figure 11B).
We observe again that the values for the fitted $Gr_{c}$ are greater than the
one proposed by Hoyal et al., (1999) by analogy with thermal convection,
whilst they also substantially differ from one another. We therefore also fit
the results to a power law $V_{f}\propto D_{p}^{\eta}$ finding $\eta=0.38\pm
0.13$ ($R^{2}=0.94$) and $\eta=0.42\pm 0.10$ ($R^{2}=0.88$) respectively to
the two volume fractions which is in very good agreement with the analytical
formulation (equation 1) that suggests $\eta=0.4$.

####  4.2.3 Particle mass flux, particle concentration in the lower layer and
accumulation rate

Given the excellent agreement between the proposed model and both LSA analysis
and analogue experiments described above, we take advantage of having 3D data
from the numerical simulations in order to extract other parameters which are
difficult to obtain otherwise (Fries et al.,, 2021). Three interesting
parameters are the particle mass flux across a plane, the particle
concentration in the lower layer and the amount of particles accumulated at
the bottom of the tank, which can be related to the accumulation rate. The
latter is especially interesting as, when the model is applied to volcanic
clouds, it could eventually be compared with field data (Bonadonna et al.,,
2011).

We calculate the mass flux across a horizontal plane (actually a thin box of
thickness $\delta x$) as shown in Figure 12A with

$J=\frac{\Delta m}{A\Delta t},$ (48)

where $\Delta m$ is the mass crossing the yellow plane of area $A$ in time
$\Delta t$, and is given by the mass difference in the volume below the plane
between $t$ and $t-\Delta t$. The mass below at each time is calculated by
summing the mass of particles in each cell $i$ of volume $\Delta V$, which is
individually given by $m_{i}=\Delta V\phi_{i}\rho_{p}$. Figure 12B shows the
temporal evolution of the particle mass flux settling through the yellow plane
(located at $0.15\,m$ below the barrier), for several initial particle volume
fractions. The vertical black dashed line indicates the theoretical time
$T_{i}$ when particles would be expected to reach the plane if they were
settling individually at their Stokes velocity. For the different simulations,
we clearly observe that the moments when the flux starts initially increasing
(i.e. the arrival of the fastest finger) are much earlier than $T_{i}$ and
this shows the extent to which the collective settling enhances the premature
sedimentation. After the initial increase, the fluxes exhibit strong
oscillations around a high plateau. These oscillations are associated with the
intermittent nature of PBL detachment and convection in the lower layer.
Finally, the particle mass flux reaches a plateau after some time which shows
the end of convection and a transition to individual settling. Throughout, the
average mass flux, as well as the amplitude of the oscillations increases with
the initial volume fraction.

Another way to highlight the enhancement of the sedimentation rate by
collective settling is to study the spatial distribution of particles beneath
the interface. Assuming a quiescent upper layer and a convective lower layer,
akin to our simulations, Hoyal et al., (1999) derived equation 4 for the
evolution of the particle concentration in the lower layer. The derivation of
this formulation assumes that $\dot{M}_{out}\neq 0$ since $t=0$ but in fact,
$\dot{M}_{out}=0$ for $t<t_{a}$ where $t_{a}$ is the time when the first
particles reach the bottom of the tank. Also, equation 4 only remains valid
for $t<t_{lim}$, where $t_{lim}=h_{1}/V_{s}$ , $h_{1}$ being the thickness of
the upper layer. After this time, there are no longer any particles remaining
in the upper layer and $\dot{M}_{in}=0$. We therefore propose an extension for
the solution of the problem (see section 2.1 in Supplementary material) which
becomes

$\begin{array}[]{rcl}C_{2}\left(t\right)=\frac{V_{s}}{h_{2}}C_{1}\left(0\right)t,&\mbox{for}&t<t_{a},\end{array}$
(49)

$\begin{array}[]{rcl}C_{2}\left(t\right)=C_{1}\left(0\right)\left[1+\left(\frac{V_{s}}{h_{2}}t_{a}-1\right)e^{-\frac{V_{s}}{h_{2}}\left(t-t_{a}\right)}\right],&\mbox{for}&t_{a}\leq
t<t_{lim},\end{array}$ (50)

$\begin{array}[]{rcl}C_{2}\left(t\right)=C_{1}\left(0\right)\left[1+\left(\frac{V_{s}}{h_{2}}t_{a}-1\right)e^{-\frac{V_{s}}{h_{2}}\left(t_{lim}-t_{a}\right)}\right]\left(1+\frac{h_{1}}{h_{2}}-\frac{V_{s}}{h_{2}}t\right),&\mbox{for}&t\geq
t_{lim},\end{array}$ (51)

where $h_{2}$ is the thickness of the lower layer. Equation 51 assumes that
the convection stops at $t_{lim}$, which suggests a quiescent settling in the
lower layer after that time with a constant flux
$\dot{M}_{out}=AV_{s}C_{2}\left(t_{l}im\right)$.

Figure 12: (A) Horizontal planar surface (yellow slice) located $0.15\,m$
below the barrier, across which the particle mass is computed in the
simulation domain. (B) Temporal evolution for the mass of particles crossing
the plane. Black dashed line: theoretical time for the particles to reach the
plane at their individual Stokes velocity.

An interesting result coming out of the previous analytical study is the mass
of particles accumulating at the bottom of the tank and the associated
accumulation rate. We can derive an analytical prediction for the mass of
particles $m_{b}$ accumulated at the bottom of the tank for the different
regimes highlighted above. Thus, by integration of the flux (see Section 2.2
in Supplementary material) we have

$\begin{array}[]{rcl}m_{b}=0,&\mbox{for}&t<t_{a},\end{array}$ (52)

$\begin{array}[]{rcl}m_{b}=m_{0}\frac{V_{s}}{h_{1}}\left[t+\left(\frac{h_{2}}{V_{s}}-t_{a}\right)e^{-\frac{V_{s}}{h_{2}}\left(t-t_{a}\right)}-\frac{h_{2}}{V_{s}}\right],&\mbox{for}&t_{a}\leq
t<t_{lim},\end{array}$ (53)

$\begin{array}[]{rcl}m_{b}=m_{0}\frac{V_{s}}{h_{1}}\\{t_{lim}-\frac{h_{2}}{V_{s}}\left[1+\left(\frac{V_{s}}{h_{2}}t_{a}-1\right)e^{-\frac{V_{s}}{h_{2}}\left(t_{lim}-t_{a}\right)}\right]\left(1+\frac{h_{1}}{h_{2}}-\frac{V_{s}}{h_{2}}t\right)\\},&\mbox{for}&t\geq
t_{lim},\end{array}$ (54)

where $m_{0}$ is the initial mass of particles injected in the upper layer.
Finally, at the time $t_{lim}+h_{2}/V_{s}$ , all the particle have settled
through the lower layer, thus $m_{b}=m_{0}$. Figure 13A shows the simulated
particle accumulation at the bottom of the tank through time, for different
particle sizes as well as the analytical prediction (equations 52-54). We
compare as well with the analytical formulation of the mass which assumes that
the lower is still turbulently convective even after the time $t_{lim}$
(equation S43 in Supplementary material, dashed lines in Figure 13A). In order
to compare between this prediction and the model results, $t_{a}$ is fitted in
order to have the best agreement between the numerical data and equations
52-54. The results show clearly that the quiescent model of the lower layer
for $t\geq t_{lim}$ agrees very well with the simulations and suggest that the
entirely convective model underestimates the accumulation rate. Additionally,
the fitted parameter $t_{a}$ is coherent with the time for the first fingers
to reach the bottom of the tank in the simulations. Figure 13B shows the
instantaneous accumulation rate computed from the numerical data for several
initial volume fractions, as estimated by

$\frac{1}{A}\frac{dm_{b}}{dt}.$ (55)

Figure 13: (A) Temporal evolution of the mass of particles accumulating at the
bottom of the tank for several particle sizes. The dashed and dotted lines
represent the extended analytical model of Hoyal et al., (1999). Particle
volume fraction $\phi_{0}=1.19\times 10^{-3}$. (B) Accumulation rate
calculated at the bottom of the tank for several particle volume fractions and
a particle size of $40\,\mu m$. The coloured dashed lines are the rate derived
from the analytical model. The black dashed line is the theoretical time at
which all particles have settled across the interface.

We observe, for each initial particle volume fraction, an initial increase of
the accumulation rate with time which reflects the enhancement of the
sedimentation process due to convection. Interestingly, the accumulation rate
then reaches a plateau at around $t=t_{lim}$, indicating that the system
switches to a steady settling regime once all particles have left the upper
layer. We compare also with the analytical relations which again have very
good agreement with our simulations. Finally, using the determined $t_{a}$, we
can also calculate the concentration $C_{2}\left(t\right)$, as calculated with
the analytical expressions in equations 49-51. Figure 14 shows a comparison
with the average $C_{2}\left(t\right)$ as measured in simulations for a
particle size of $40\,\mu m$ and three different initial upper layer
concentrations, finding very good agreement.

Figure 14: Evolution of the average particle volume fraction in the lower
layer for particle of size $40\,\mu m$. Black: $C_{0}=2g/l$
$(\phi_{0}=7.94\times 10^{-4})$, Red: $C_{0}=4g/l$ $(\phi_{0}=1.59\times
10^{-3})$ and Blue: $C_{0}=6g/l$ $(\phi_{0}=2.38\times 10^{-3})$. Solid lines:
numerical model. Dashed lines: modified Hoyal et al., (1999) model (equations
4 and 49-51).

## 5\. Discussion

### 5.1. Model caveats

Our numerical model has been validated by comparing various outputs with
results from linear stability analysis, lab experiments (Fries et al.,, 2021)
and theoretical predictions from previous studies (Carazzo and Jellinek,,
2012; Hoyal et al.,, 1999). Even though these comparisons are good (Figures
9-14), the results provided by the model inherits the caveats of the
experiments. Indeed, the static and confined configuration, as well as the
fact that we performed the simulations in water, mean that we cannot fully
extend the results to the volcanic case yet. Thus, further investigations are
necessary to better simulate the volcanic environment (in air, presence of
wind…). Additionally, it is necessary to consider the limits of validity of
the assumption that the particles can be represented as a continuum. Whilst
the condition on the particle coupling is given by the Stokes number ($St<1$),
there is also a condition on the particle volume fraction to take in account.
Harada et al., (2013) and Yamamoto et al., (2015) derived a dimensionless
number in order to characterise the transition between fluid-like and
particle-like settling. Although this number is only valid for narrow channel
configurations, which are considerably different from ours, it highlights the
fact that the particle size, volume fraction and characteristic length scale
of the flow are critical parameters to define the validity of the continuum
assumption. Thus, the transition from fluid-like to particle-like behaviour is
achieved by decreasing the volume fraction and characteristic length scale and
increasing the particle size. Near this transition, the use of a single-phase
model, such as that presented here, should be treated with caution and this
reveals the need for a comparison with future models which explicitly account
for the drag contribution of individual particles.

Another related caveat concerns the numerical diffusion underlying the use of
a Eulerian approach to describe the transport of particles. Compared to
classical first order finite difference methods, the use of the third order
WENO procedure has drastically reduced the numerical diffusion. It is also
possible to further reduce the induced numerical diffusion by increasing the
order of the WENO scheme (i.e. increase also the computational cost). However,
for problems purely related to advection, where the presence of any diffusion
is critical, another strategy, such as two-phase models (using a Lagrangian
approach where individual particles are explicitly modelled), has to be
considered.

### 5.2. Vertical finger velocity

We have compared the simulated vertical velocity of fingers with experimental
observations (Fries et al.,, 2021) and a theoretical prediction (equation 1)
from (Carazzo and Jellinek,, 2012; Hoyal et al.,, 1999) (Figure 11). This
expression depends on a critical Grashof number, which by analogy with thermal
convection (Turner,, 1973) has previously been assumed to be $10^{3}$ (Hoyal
et al.,, 1999). This value effectively corresponds to a dimensionless critical
PBL thickness at which point the PBL can detach and form fingers. However,
both the model results and experimental observations summarised in Figure 11
suggest that $Gr_{c}>10^{3}$ for our configuration. Furthermore, as seen in
Figure 11B, the curve for $V_{f}$ using $Gr_{c}=10^{3}$ (blue dotted line)
crosses the Stokes velocity curve around $95\,\mu m$ for instance with an
initial particle volume fraction of $\phi_{0}=1.19\times 10^{-3}$, suggesting
this value should be the upper particle size limit for finger formation.
However, in agreement with experiments (Fries et al.,, 2021; Scollo et al.,,
2017), we observe a larger threshold for the finger formation to be in the
size range $[115-125]\,\mu m$, for $\phi_{0}=1.19\times 10^{-3}$, and in the
range $[145-160]\,\mu m$ for $\phi_{0}=3.57\times 10^{-3}$, in this particular
configuration. We also showed that equation 1 poorly predicts the observed
dependence of the finger velocity on the initial particle volume fraction.
Indeed, our studies suggests an alternative power law that better describes
the dependence of $V_{f}$ on $\phi_{0}$. Equation 1 has been derived by a
scaling theory that involves $\delta_{PBL}$ as characteristic length of the
problem (Carazzo and Jellinek,, 2012; Hoyal et al.,, 1999) and the
discrepancies highlighted in this paper (Figure 11) may suggest that
$\delta_{PBL}$ actually has a slightly different dependence on the initial
particle volume fraction. Moreover, the use of the Grashof number as an
appropriate scaling for the PBL thickness remains uncertain. On one hand, our
results suggest that if instability does occur once a critical Grashof number
is reached, the critical value taken from the thermal convection analogy is
not valid. On the other hand, the Grashof number may simply not be the correct
dimensionless form of the PBL thickness, and different flow configurations
will produce different critical values. However, the predicted dependence of
the finger velocity on the particle diameter by equation 1 shows a very good
agreement with our simulated results, as confirmed by a power-law fitting
between $V_{f}$ and $D_{p}$. Thus, whilst we have demonstrated the need for a
better scaling of $\delta_{PBL}$, equation 1 can still provide a good estimate
for the particle size threshold to form fingers. Consequently, if the size
threshold to form fingers is given when equation 1 equals the Stokes velocity
(equation 13) we can derive a formulation for the threshold

$D^{*}_{p}=\left[\frac{\left(18\mu\right)^{2}\phi\delta_{PBL}}{g\left(\rho_{p}-\rho_{f}\right)\rho_{f}}\sqrt{\left(\frac{\pi}{4}\right)}\right]^{\frac{1}{4}},$
(56)

The main caveats for this formulation are that it strongly depends on having a
correct scaling for $\delta_{PBL}$ and obviously, this estimation is valid
under the assumption that particles settle at their Stokes velocity, which is
reasonable for our study but might be uncertain in nature where the ambient
fluid is air and for non-spherical particles.

### 5.3. Particle concentration in the lower layer and mass accumulation rate

We have proposed a modified analytical formulation for the particle
concentration in the lower layer $C_{2}\left(t\right)$ and consequently for
the mass of particles accumulated at the bottom of the tank
$m_{b}\left(t\right)$. Despite some numerical artefacts that can be seen on
Figure 13B where the computed accumulation rate seems to be non-zero before
$t_{a}$, there is very good agreement between the simulations and the
analytical model. The artefacts themselves are due to fluctuating numerical
errors that do not affect the final results.

The analytical predictions for $C_{2}\left(t\right)$ and $m_{b}\left(t\right)$
are step-wise functions depending on $t_{a}$, the time it takes for the first
particles to reach the bottom. For $t<t_{a}$, the analytical model predicts
that $C_{2}\left(t\right)$ increases linearly with time since the formulation
assumes that, during this period, particles are settling individually. In
fact, our numerical results show that convective settling does occur for
$t<t_{a}$ but, since this time period is short, the linear law seems to be a
satisfactory approximation for the early-time average lower layer particle
concentration. However, in order to compare our simulated results with the
analytical prediction, we fitted the parameter $t_{a}$ in this study. Although
we are able to obtain excellent agreement between model and theory, it would
be better to develop a fully independent formulation. To achieve this, it is
necessary to also provide an analytical estimation for $t_{a}$. One possible
approach would be to assume the decomposition
$t_{a}=t_{a}^{{}^{\prime}}+t_{a}^{{}^{\prime\prime}}$ where
$t_{a}^{{}^{\prime}}$ is the time during which the PBL initially grows beneath
the interface at the individual particle settling velocity, i.e.,
$t_{a}^{{}^{\prime}}=\delta_{PBL}/V_{s}$ , and $t_{a}^{{}^{\prime\prime}}$ is
the time between the PBL detachment and the first arrival of particles at the
base of the domain. If, during this stage, we assume that the particles are
advected at the finger velocity then
$t_{a}^{{}^{\prime\prime}}=\left(h_{2}-\delta_{PBL}\right)/V_{f}$. We
therefore see that $t_{a}$ strongly depends on $\delta_{PBL}$, which
highlights once again the need for a correct scaling of the PBL thickness, as
discussed in the previous section. Another interesting result concerns the
accumulation rate of particles at the base of the domain in the presence of
fingers. Figure 13B shows the accumulation rate increases with time for
$t_{a}<t<t_{lim}$, in agreement with the analytical prediction (i.e.
combination of equations 53 and 55 which provides an exponential increase of
$m_{b}$). Conversely, if the particles had settled individually, the
accumulation rate would be temporally constant. This shows that temporally
resolved measurements of the accumulation rate of particles from volcanic
clouds may record temporal signatures of sedimentation via settling-driven
gravitational instabilities. Whilst there is already a spatial deposit
signature of settling-driven gravitational instabilities (i.e. bimodal
grainsize distribution) (Bonadonna et al.,, 2011; Manzella et al.,, 2015),
this is not unique and can be generated by other mechanisms such as particle
aggregation (Brown et al.,, 2012). Accumulation rate data from the field may
therefore provide a powerful tool for distinguishing the efficiency of
convective sedimentation beneath volcanic clouds.

## 6\. Conclusions

We have presented an innovative hybrid Lattice Boltzmann Finite Difference 3D
model in order to simulate settling-driven gravitational instabilities at the
base of volcanic ash clouds. Such instabilities occur when particles settle
through a density interface at the base of a suspension, leading to the
formation of an unstable particle boundary layer (Carazzo and Jellinek,, 2012;
Hoyal et al.,, 1999; Manzella et al.,, 2015), and also occur in other natural
settings, such as river plumes (Davarpanah Jazi and Wells,, 2016). Our
numerical model makes use of a low-diffusive WENO procedure to solve the
advection-diffusion-settling equation for the particle volume fraction. The
use of such a routine allows us to minimise errors associated with numerical
diffusion and has the advantage of being applied to simple uniform meshes,
which makes the coupling with the LBM easier. This innovative use of the WENO
scheme, therefore, represents an effective tool for the solving of advection-
dominated problems. Our implementation of the third order WENO finite
difference scheme will be integrated in a future release of the open-source
$Palabos$ code. Our model has been successfully validated by comparing the
results with i) predictions from linear stability analysis where we show that
the model is able to simulate settling-driven gravitational instabilities from
the initial disturbance through the linearly-unstable regime, ii) analogue
experiments (Fries et al.,, 2021) and iii) theoretical models (Carazzo and
Jellinek,, 2012; Hoyal et al.,, 1999) in order to reproduce the non-linear
regime which describes the downward propagation of fingers. We also confirmed
the premature sedimentation process through collective settling compared to
individual settling.

Our model provides new insights into:

  * •

the value of the critical Grashof number. From measurements of the vertical
finger speed, we have found $Gr_{c}\sim 10^{4}$ in our configuration. This
value differs from the one suggested by analogy with thermal convection
($Gr_{c}\sim 10^{3}$) (Hoyal et al.,, 1999). Our results suggest that either
the critical Grashof number for settling-driven gravitational instabilities is
greater than in the thermal convection case or that the Grashof number may not
be the correct dimensionless form of the PBL thickness. In any case, this
highlights the need for further investigation of the scaling of the PBL
thickness $\delta_{PBL}$.

  * •

the presence of a particle size threshold for the finger formation. Using our
results, we have proposed an analytical formulation for this threshold
depending on the density of particles, the viscosity of the medium and also
the bulk density difference between the two fluid layers.

  * •

the signature of settling-driven gravitational instabilities (i.e.
accumulation rate). We show that the accumulation rate of particles at the
tank base initially increases with time before reaching a plateau. This
contrasts with the constant accumulation rate associated with individual
particle settling. This suggests that accumulation rate data could be used
during tephra fallout to distinguish between sedimentation through settling-
driven gravitational instabilities and individual-particle sedimentation.

We have also demonstrated how our numerical model can be used to expand the
initial conditions and configuration settings that can be explored through
experimental investigations. The results presented so far in an aqueous media
permitted model validation but have also opened fundamental questions that
will be addressed in future works involving configurations more similar to the
natural system. Indeed, thanks to the strengths of the LBM, the model can
easily be applied to more complex systems and provide a robust tool for the
transition from the laboratory studies to volcanic systems, as well as other
environmental flows.

## Conflict of Interest

The authors declare that the research was conducted in the absence of any
commercial or financial relationships that could be construed as a potential
conflict of interest.

## Data Availability Statement

The datasets generated for this study can be found in the
article/Supplementary Material. Additional raw data are available upon
request.

## Author Contributions

Jonathan Lemus integrated the WENO procedure in the Palabos framework,
conducted the simulations and data analysis under the supervision of Paul
Jarvis, Jonas Lätt, Costanza Bonadonna and Bastien Chopard. Jonathan Lemus
drafted the manuscript. Jonas Lätt and Bastien Chopard were involved in the
development of the Palabos code. All authors have contributed to data
interpretation as well as the editing and finalising of the paper.

## Funding

The study has been funded by the Swiss National Science Foundation $\sharp
200021\\_169463$.

## Acknowledgments

All the simulations presented in this paper have been performed using the High
Performance Computing (HPC) facilities $Baobab$ and $Yggdrasil$ of the
University of Geneva. We would like to thank Amanda B. Clarke and Jeremy C.
Phillips for constructive discussions about the problem.

## References

  * Alsinan et al., (2017) Alsinan, A., Meiburg, E., and Garaud, P. (2017).  A settling-driven instability in two-component, stably stratified fluids.  Journal of Fluid Mechanics, 816:243–267.
  * Baltensperger and Berrut, (2001) Baltensperger, R. and Berrut, J.-p. (2001).  The linear rational collocation method.  134:243–258.
  * Berrut and Mittelmann, (2004) Berrut, J. P. and Mittelmann, H. D. (2004).  Adaptive point shifts in rational approximation with optimized denominator.  Journal of Computational and Applied Mathematics, 164-165:81–92.
  * Bhatnagar et al., (1954) Bhatnagar, P., Gross, E., and Krook, M. (1954).  A model for collision processes in gases.  Physical Review, 94(3):515–525.
  * Blong, (2000) Blong, R. (2000).  Volcanic hazards and risk management.  In Encyclopedia of volcanoes, pages 1215–1227. Academic Press.
  * Bonadonna et al., (2021) Bonadonna, C., Biass, S., Menoni, S., and Gregg, C. E. (2021).  Assessment of risk associated with tephra-related hazards.
  * Bonadonna et al., (2012) Bonadonna, C., Folch, A., Loughlin, S., and Puempel, H. (2012).  Future developments in modelling and monitoring of volcanic ash clouds: Outcomes from the first IAVCEI-WMO workshop on Ash Dispersal Forecast and Civil Aviation.  Bulletin of Volcanology, 74(1):1–10.
  * Bonadonna et al., (2011) Bonadonna, C., Genco, R., Gouhier, M., Pistolesi, M., Cioni, R., Alfano, F., Hoskuldsson, A., and Ripepe, M. (2011).  Tephra sedimentation during the 2010 Eyjafjallajkull eruption (Iceland) from deposit, radar, and satellite observations.  Journal of Geophysical Research: Solid Earth, 116(12).
  * Brown et al., (2012) Brown, R. J., Bonadonna, C., and Durant, A. J. (2012).  A review of volcanic ash aggregation.  Physics and Chemistry of the Earth, 45-46:65–78.
  * Burgisser et al., (2005) Burgisser, A., W. Bergantz, G., and E. Breidenthal, R. (2005).  Addressing Complexity in Laboratory Experiments: The Scaling of Dilute Multiphase Flows in Magmatic Systems.  Journal of Volcanology and Geothermal Research, 141:N{$\backslash$textdegree} 3–4, 245–265.
  * Burns and Meiburg, (2012) Burns, P. and Meiburg, E. (2012).  Sediment-laden fresh water above salt water: Linear stability analysis.  Journal of Fluid Mechanics, 691:279–314.
  * Burns and Meiburg, (2014) Burns, P. and Meiburg, E. (2014).  Sediment-laden fresh water above salt water: Nonlinear simulations.  Journal of Fluid Mechanics, 762:156–195.
  * Carazzo and Jellinek, (2012) Carazzo, G. and Jellinek, A. M. (2012).  A new view of the dynamics, stability and longevity of volcanic clouds.  Earth and Planetary Science Letters, 325-326:39–51.
  * Cardoso and Zarrebini, (2001) Cardoso, S. S. S. and Zarrebini, M. (2001).  Convection driven by particle settling surrounding a turbulent plume.  Chemical Engineering Science, 56(11):3365–3375.
  * Carey, (1997) Carey, S. (1997).  Influence of convective sedimentation on the formation of widespread tephra fall layers in the deep sea.  Geology, 25(9):839–842.
  * Chandrasekhar, (1961) Chandrasekhar, S. (1961).  Hydrodynamic and hydromagnetic stability.  Oxford edition.
  * Chen and Doolen, (2010) Chen, S. and Doolen, G. D. (2010).  Lattice Boltzmann Method for fluid flows.  Scholarpedia, 5(5):9507.
  * Chou and Shao, (2016) Chou, Y. J. and Shao, Y. C. (2016).  Numerical study of particle-induced Rayleigh-Taylor instability: Effects of particle settling and entrainment.  Physics of Fluids, 28(4).
  * Crimaldi, (2008) Crimaldi, J. P. (2008).  Planar laser induced fluorescence in aqueous flows.  Experiments in Fluids, 44(6):851–863.
  * Davarpanah Jazi and Wells, (2016) Davarpanah Jazi, S. and Wells, M. G. (2016).  Enhanced sedimentation beneath particle-laden flows in lakes and the ocean due to double-diffusive convection.  Geophysical Research Letters, 43(20):10,883–10,890.
  * Durant, (2015) Durant, A. J. (2015).  RESEARCH FOCUS: Toward a realistic formulation of fine-ash lifetime in volcanic clouds.  Geology, 43:271–272.
  * Folch, (2012) Folch, A. (2012).  A review of tephra transport and dispersal models : Evolution , current status , and future perspectives.  Journal of Volcanology and Geothermal Research, 235-236:96–115.
  * Folch et al., (2020) Folch, A., Mingari, L., Gutierrez, N., Hanzich, M., Macedonio, G., and Costa, A. (2020).  FALL3D-8.0: A computational model for atmospheric transport and deposition of particles, aerosols and radionuclides - Part 1: Model physics and numerics.  Geoscientific Model Development, 13(3):1431–1458.
  * Fries et al., (2021) Fries, A., Lemus, J., Jarvis, P. A., Clarke, A. B., Phillips, J. C., Manzella, I., and Bonadonna, C. (2021).  The Influence of Particle Concentration on the Formation of Settling-Driven Gravitational Instabilities at the Base of Volcanic Clouds.  Frontiers Earth Sciences.
  * Godunov, (1954) Godunov (1954).  Different methods for shock waves.  Moscow State University.
  * Godunov, (1959) Godunov (1959).  A difference scheme for numerical solution of discontinuous solution of hydrodynamic equations.  Math. Sbornik, 47:271–306.
  * Gouhier et al., (2019) Gouhier, M., Eychenne, J., Azzaoui, N., Guillin, A., Deslandes, M., Poret, M., Costa, A., and Husson, P. (2019).  Low efficiency of large volcanic eruptions in transporting very fine ash into the atmosphere.  Scientific Reports, 9(1):1–12.
  * Gudmundsson, (2011) Gudmundsson, G. (2011).  Respiratory health effects of volcanic ash with special reference to Iceland. A review.  Clinical Respiratory Journal, 5(1):2–9.
  * Guffanti et al., (2008) Guffanti, M., Mayberry, G., Casadevall, T., and Wunderman, R. (2008).  Volcanic hazards to airports.  Natural Hazards, 51:287–302.
  * Guo et al., (2002) Guo, Z., Zheng, C., and Shi, B. (2002).  Discrete lattice effects on the forcing term in the lattice Boltzmann method.  65:1–6.
  * Harada et al., (2013) Harada, S., Kondo, M., Watanabe, K., Shiotani, T., and Sato, K. (2013).  Collective settling of fine particles in a narrow channel with arbitrary cross-section.  Chemical Engineering Science, 93:307–312.
  * He and Luo, (1997) He, X. and Luo, L.-s. (1997).  Lattice Boltzmann Model for the Incompressible Navier-Stokes Equation.  88:927–944.
  * Hosseini et al., (2017) Hosseini, S. A., Darabiha, N., Thévenin, D., and Eshghinejadfard, A. (2017).  Stability limits of the single relaxation-time advection-diffusion lattice Boltzmann scheme.  International Journal of Modern Physics C, 28(12):1–31.
  * Hoyal et al., (1999) Hoyal, D., Bursik, I., and Atkinson, F. (1999).  Settling-driven convection: A mechanism of sedimentation from stratified fluids.  104:7953–7966.
  * Jacobs et al., (2013) Jacobs, C. T., Collins, G. S., Piggott, M. D., Kramer, S. C., and Wilson, C. R. (2013).  Multiphase flow modelling of volcanic ash particle settling in water using adaptive unstructured meshes.  Geophysical Journal International, 192(2):647–665.
  * Jenkins et al., (2015) Jenkins, S. F., Wilson, T., Magill, C., Miller, V., Stewart, C., Blong, R., Marzocchi, W., Boulton, M., Bonadonna, C., and Costa, A. (2015).  Volcanic ash fall hazard and risk.  Number 2015.
  * Jiang and Shu, (1996) Jiang, G.-S. and Shu, C.-W. (1996).  Efficient Implementation of Weighted ENO Schemes.  Journal of Computational Physics, 228(126):202–228.
  * Jones et al., (2007) Jones, A., Thomson, D., Hort, M., and Devenish, B. (2007).  The U.K. Met Office’s Next-Generation Atmospheric Dispersion Model, NAME III.  Air Pollution Modeling and Its Application XVII, (December 2015):580–589.
  * Keck et al., (2021) Keck, J. B., Cottet, G. H., Meiburg, E., Mortazavi, I., and Picard, C. (2021).  Double-diffusive sedimentation at high Schmidt numbers: Semi-Lagrangian simulations.  Physical Review Fluids, 6(2):1–10.
  * Koochesfahani, (1984) Koochesfahani, M. M. (1984).  Experiments on turbulent mixing and chemical reactions in a liquid mixing layer.  PhD thesis, California Institute of Technology.
  * Kruger et al., (2017) Kruger, T., Kusumaatmaja, H., Kuzmin, A., Shardt, O., Goncalo, S., and Viggen, E. M. (2017).  The lattice boltzmann method, principles and practice.  Number 207.
  * Latt et al., (2020) Latt, J., Malaspinas, O., Kontaxakis, D., Parmigiani, A., Lagrava, D., Brogi, F., Belgacem, M. B., Thorimbert, Y., Leclaire, S., Li, S., Marson, F., Lemus, J., Kotsalos, C., Conradin, R., Coreixas, C., Petkantchin, R., Raynaud, F., Beny, J., and Chopard, B. (2020).  Palabos: Parallel Lattice Boltzmann Solver.  Computers and Mathematics with Applications, (xxxx).
  * Leclaire et al., (2017) Leclaire, S., Parmigiani, A., Malaspinas, O., Chopard, B., and Latt, J. (2017).  Generalized three-dimensional lattice Boltzmann color-gradient method for immiscible two-phase pore-scale imbibition and drainage in porous media.  Physical Review E, 95(3):33306.
  * Lewis, (1950) Lewis, D. J. (1950).  The instability of liquid surfaces when accelerated in a direction perpendicular to their planes. II.  Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 202(1068):81–96.
  * Liu et al., (1994) Liu, X.-D., Osher, S., and Chan, T. (1994).  Weighted Essentially Non-Oscillatory Schemes.  Journal of Computational Physics, pages 1–32.
  * Manville and Wilson, (2004) Manville, V. and Wilson, C. J. (2004).  Vertical density currents: A review of their potential role in the deposition and interpretation of deep-sea ash layers.  Journal of the Geological Society, 161(6):947–958.
  * Manzella et al., (2015) Manzella, I., Bonadonna, C., Phillips, J. C., and Monnard, H. (2015).  The role of gravitational instabilities in deposition of volcanic ash.  Geology, 43(3):211–214.
  * Noriega et al., (2013) Noriega, H., Reggio, M., and Vasseur, P. (2013).  Natural convection of nanofluids in a square cavity heated from below.  Computational Thermal Sciences: An International Journal, vol. 5, no.
  * Parmigiani et al., (2009) Parmigiani, A., Huber, C., Chopard, B., Latt, J., and Bachmann, O. (2009).  Application of the multi distribution function lattice Boltzmann approach to thermal flows.  The European Physical Journal Special Topics, 171(1):37–43.
  * Prata et al., (2021) Prata, A. T., Mingari, L., Folch, A., Macedonio, G., and Costa, A. (2021).  FALL3D-8.0: A computational model for atmospheric transport and deposition of particles, aerosols and radionuclides - Part 2: Model validation.  Geoscientific Model Development, 14(1):409–436.
  * Prata and Rose, (2015) Prata, F. and Rose, B. (2015).  Volcanic Ash Hazards to Aviation.  pages 911–934. Academic Press, Amsterdam.
  * Roche and Carazzo, (2019) Roche, O. and Carazzo, G. (2019).  The contribution of experimental volcanology to the study of the physics of eruptive processes, and related scaling issues: A review.  Journal of Volcanology and Geothermal Research, 384:103–150.
  * Scollo et al., (2017) Scollo, S., Bonadonna, C., and Manzella, I. (2017).  Settling-driven gravitational instabilities associated with volcanic clouds: new insights from experimental investigations.  Bulletin of Volcanology, 79(6):1–14.
  * Scollo et al., (2010) Scollo, S., Folch, A., Coltelli, M., and Realmuto, V. J. (2010).  Three-dimensional volcanic aerosol dispersal: A comparison between Multiangle Imaging Spectroradiometer (MISR) data and numerical simulations.  Journal of Geophysical Research Atmospheres, 115(24):1–14.
  * Scollo et al., (2008) Scollo, S., Tarantola, S., Bonadonna, C., Coltelli, M., and Saltelli, A. (2008).  Sensitivity analysis and uncertainty estimation for tephra dispersal models.  Journal of Geophysical Research: Solid Earth, 113(6):1–17.
  * Sharp, (1984) Sharp, D. H. (1984).  An overview of Rayleigh-Taylor instability.  Physica D: Nonlinear Phenomena, 12(1-3):3–18.
  * Spence et al., (2005) Spence, R. J. S., Kelman, I., Baxter, P. J., Zuccaro, G., and Petrazzuoli, S. (2005).  Natural Hazards and Earth System Sciences Residential building and occupant vulnerability to tephra fall.  Natural Hazards and Earth System Sciences, 5:477–494.
  * Stokes, (1851) Stokes, G. G. (1851).  On the Theories of the Internal Friction of Fluids in Motion, and of the Equilibrium and Motion of Elastic Solids.  Mathematical and Physical Papers vol.1, 9:75–129.
  * Turner, (1973) Turner, J. S. (1973).  Buoyancy Effects in Fluids.  Cambridge University Press, Cambridge.
  * Völtz et al., (2001) Völtz, C., Pesch, W., and Rehberg, I. (2001).  Rayleigh-Taylor instability in a sedimenting suspension.  Physical Review E, 65(1):011404.
  * Watt et al., (2015) Watt, S. F., Gilbert, J. S., Folch, A., Phillips, J. C., and Cai, X. M. (2015).  An example of enhanced tephra deposition driven by topographically induced atmospheric turbulence.  Bulletin of Volcanology, 77(5).
  * Wilson et al., (2011) Wilson, T. M., Cole, J. W., Stewart, C., Cronin, S. J., and Johnston, D. M. (2011).  Ash storms: Impacts of wind-remobilised volcanic ash on rural communities and agriculture following the 1991 Hudson eruption, southern Patagonia, Chile.  Bulletin of Volcanology, 73(3):223–239.
  * Yamamoto et al., (2015) Yamamoto, Y., Hisataka, F., and Harada, S. (2015).  Numerical simulation of concentration interface in stratified suspension: Continuum-particle transition.  International Journal of Multiphase Flow, 73:71–79.
  * Yu et al., (2013) Yu, X., Hsu, T. J., and Balachandar, S. (2013).  Convective instability in sedimentation: Linear stability analysis.  Journal of Geophysical Research: Oceans, 118(1):256–272.
  * Yu et al., (2014) Yu, X., Hsu, T.-J., and Balachandar, S. (2014).  Convective instability in sedimentation: 3-D numerical study.  Journal of Geophysical Research: Oceans, (119):3868–3882.