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

Field observations and laboratory experiments have shown that ash sedimentation can be significantly affected by collective settling mechanisms that promote premature ash deposition, with important implications for associated impacts. Among these mechanisms, settling-driven gravitational instabilities result from the formation of a gravitationally-unstable particle boundary layer (PBL) that grows between volcanic ash clouds and the underlying atmosphere. The PBL destabilises once it reaches a critical thickness, triggering the formation of rapid, downward-moving ash fingers that remain poorly characterised. We simulate this process by coupling a Lattice Boltzmann model, which solves the Navier-Stokes equations for the fluid phase, with a Weighted Essentially Non Oscillatory (WENO) finite difference scheme which solves the advection-diffusion-settling equation describing particle transport. Since the physical problem is advection dominated, the use of the WENO scheme reduces numerical diffusivity and ensures accurate tracking of the temporal evolution of the interface between the layers. We have validated the new model by showing that the simulated early-time growth rate of the instability is in very good agreement with that predicted by linear stability analysis, whilst the modelled late-stage behaviour also successfully reproduces quantitative results from published laboratory experiments.


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 Modelling gravitational instabilities 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 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).  (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 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;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 , 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;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; where ρ PBL is the PBL bulk density, ρ f the underlying fluid density, g = 9.81m.s −2 the gravitational acceleration and δ PBL the PBL thickness, which by analogy with thermal convection (Turner, 1973) is taken to be  δ PBL = Gr c ν 2 g where g = g ρ PBL − ρ f /ρ f , ν 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 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 , although recent experimental observations suggests Gr c ≈ 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 supV s . 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,  also developed a series of analytical massbalance 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Ṁ in and the mass flux of particle leaving by sedimentationṀ out where t is time. Assuming that M 2 (t) = Ah 2 C 2 (t), where C 2 (t) is the average particle concentration in the lower layer,  solved this equation usingṀ in = AV s C 1 (0),Ṁ out = AV s C 2 (t) and the initial condition C 2 (0) = 0. Thus where C 1 (0) 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) where ρ p is the particle density, D p the particle diameter, µ 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.

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 φ ( x,t) (Burns and Meiburg, 2014;Chou and Shao, 2016;Yu et al., 2014). φ ( x,t) satisfies the advection-diffusion-settling equation where u f ( x,t) is the fluid velocity, D c the particle diffusion coefficient, e z the vertical unit vector and x = (x, y, z) the position coordinate. The fluid is considered incompressible meaning ∇. u f = 0. Thus, equation 6 becomes The particle settling velocity depends on the ambient fluid density ρ, 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 where ρ 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 where p is the pressure and F the buoyant body force term. We complete the system of equations by taking this force term to be a function of φ and ρ 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 φ i , where each one is associated with a different size class and individually satisfies equation 7. Furthermore, the body force term becomes

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 ± 0.5 µm (measured using laser diffraction with a Bettersizer S3 Plus) and a density ρ p of 2519.4 ± 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 (0) of the upper layer is varied from 1 to 10 g/l (see Table 1 for the conversion to particle volume fraction φ 0 ). The lower layer density is kept constant at ρ 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.

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)  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 ρ 0 where S is the sugar concentration and ρ = ρ 0 (1 + αS), with α 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 φ and ρ to avoid any flux of particles or sugar across the walls. Thus we impose and on the wall nodes. Furthermore, we define the following initial states for φ and S and where φ 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 u f ( x,t = 0) = 0.

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 φ and S.

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 where δt is the time step, f eq i ρ, u f the equilibrium distribution function, τ the relaxation time associated with the flow viscosity and 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. 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 and whilst the kinematic viscosity controls the relaxation to equilibrium through the relaxation time The variable c 2 s is commonly called the speed of sound and is equal to (1/3)(δ x 2 /δt 2 ) where δ 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 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 The determination of the coefficient b, as well as the power series expansion of F i are described by . 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.

Transport of particles and other density-altering quantities
The particles and other density-altering quantities are described by continuum fields that follow an advectiondiffusion 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 ν However, a stability condition for a LBM-BGK algorithm is τ/δt = 1/2. Thus, since the problem is convection dominated, the low diffusion coefficient (D 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 .
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 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 φ , which includes an additional velocity contribution due to the settling, the truncation error associated with the first-order scheme becomes non-negligible. Whilst decreasing δ 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 . Therefore, we choose here to implement the low diffusive WENO procedure for the solution of φ , 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.

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.

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.

Nondimensionalisation
We nondimensionalise our system of equations by defining and where l c , t c and p c are characteristic quantities. We also define the dimensionless parameters and noting that Fr is a Froude number and Sc i are Schmidt numbers. Furthermore, the stream function ψ is defined such that u f = (∂ ψ/∂ z, −∂ ψ/∂ x) and the vorticity as ω = ∇ × 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) and Note that here we have neglected the term −φ ∇. (V s e z ) in equation 7 assuming the fluid density variation across the interface is sufficiently small that it does not affect the particle settling velocity.

Variable expansion and eigenvalue problem
We linearise the system of equations by expanding each variable in terms of a base state and a perturbation where ϕ (y, z,t) = {ψ, ω, φ , S},φ (z) = {ψ,ω,φ ,S} the associated base state and ϕ (y, z,t) = {ψ , ω , φ , S } the perturbation. We choose the following base states where z φ (T ) and z S (T ) 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 whereφ (z) is the perturbation amplitude, k the wavenumber and σ the instability growth rate. The linearised system of equations is then formulated in matrix form so that the problem is reduced to the eigenvalue and, in a reference frame moving downward at V s , the matrices K and W are given by and where D z = ∂ /∂ z, M = −k 2 + D 2 z and I is the identity operator. The eigenvalues σ 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 σ 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.

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;.

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 (y, z) 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 φ = φ 0 /2 and also define H (y) 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., φ 0 /10 or φ 0 /5 (relative change 3%). Figure  5A shows an example of a space-time diagram showing the evolution of H (y) through time. Furthermore, Figure 4: Dispersion relation obtained from LSA for several initial particle concentrations.
by calculating the Fourier transformH (k,t) of H (y) 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) Γ H (k,t) = (1/ (k S L S )) |H (k,t) | 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 Γ H (k,t) and compare the results with the predicted growth rates from LSA. We apply this analysis during a period when the amplitude |H (k,t) | 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) with |H i (k) | the initial amplitude and σ sim (k) the instability growth rate as determined from the simulations. Thus, the PSD can be expressed as where Γ H i is the initial spectral density. 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 (φ 0 = 1.59 × 10 −3 , 2.38 × 10 −3 and 3.97 × 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 Γ H (k sim ,t), we apply our exponential fitting (equation 46) to determine the growth rate σ 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 σ 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 φ (T ) and z S (T ) (equations 39-40) and thus determine the base states ofφ (z, T ) andS(z, T ) (Figure 7). We then perform the LSA for each φ 0 , using the appropriate base states, and obtain a dispersion relation σ = f (k). Using this relation, we predict the different growth rates σ = σ LSA,i associated with k = k sim,i and we compare with σ (k sim,i ) as measured in our simulations. Figure 8 shows the comparison between σ sim (black dots) and σ 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 σ LSA,i on the initial particle concentration is in good agreement with the simulated growth rate.   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.

Characterisation of the PBL and effect of the initial particle volume fraction on the finger velocity
The bulk density profile ρ blk , derived from the contributions of the particle concentration and sugar profiles, is given everywhere by the relation Figure 10 shows the profiles of φ , ρ f and ρ blk in the numerical simulations as well as in the experiments 8 seconds after the barrier removal for the same initial conditions (φ 0 = 3.18 × 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 φ 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 φ 0 is non-linear. By analogy with thermal convection, it has previously been assumed that Gr c = 10 3 , 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 ± 0.4 × 10 4 (R 2 = 0.92), which is an order of magnitude higher than the value previously assumed (Carazzo and Jellinek, 2012;. This agrees reasonably with the experiments, where the best fit is obtained for Gr c = 1.9 ± 0.7 × 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 φ 0 , assuming a power law of the form V f ∝ φ q 0 . According to equation 1, q = 4/15 ≈ 0.27. However, from the experiments, we obtain q = 0.50 ± 0.16 (with R 2 = 0.95) while for our simulations q = 0.37 ± 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.

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 µm (for φ 0 = 1.19×10 −3 ) and 145 µm (for φ 0 = 3.57×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 µm and 145 µm respectively for φ 0 = 1.19 × 10 −3 and φ 0 = 3.57 × 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 ∼ 125 µm with in initial particle volume fraction of φ 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 ± 3.6 × 10 3 (with R 2 = 0.91) and Gr c = 2.7 ± 0.8 × 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  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 ∝ D η p finding η = 0.38 ± 0.13 (R 2 = 0.94) and η = 0.42 ± 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 η = 0.4.

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 δ x) as shown in Figure 12A with where ∆m is the mass crossing the yellow plane of area A in time ∆t, and is given by the mass difference in the volume below the plane between t and t − ∆t. The mass below at each time is calculated by summing the mass of particles in each cell i of volume ∆V , which is individually given by m i = ∆V φ i ρ 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,  derived equation 4 for the evolution of the particle concentration in the lower layer. The derivation of this formulation assumes thatṀ out = 0 since t = 0 but in fact,Ṁ 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Ṁ in = 0. We therefore propose an extension for the solution of the problem (see section 2.1 in Supplementary material) which becomes 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Ṁ out = AV s C 2 (t l im). whilst for the experiments Gr c = 1.9 ± 0.7 × 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 φ 0 = 1.19 × 10 −3 , the best fit gives Gr c = 7.6 × 10 3 and no fingers are observed to form for particle sizes higher than 115 µm. For φ 0 = 3.57 × 10 −3 , the best fit gives Gr c = 2.7 × 10 4 and no fingers are observed to form for particle sizes higher than 145 µm. In the two plots, the blue dashed line shows equation 1 using Gr c = 10 3 from the analogy with thermal convection ).
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 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 ≥ 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 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 (t), as calculated with the analytical expressions in equations 49-51. Figure 14 shows a comparison with the average C 2 (t) as measured in simulations for a particle size of 40 µm and three different initial upper layer concentrations, finding very good agreement.

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;. 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.

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; (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 . 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 µm for instance with an initial particle volume fraction of φ 0 = 1.19 × 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] µm, for φ 0 = 1.19 × 10 −3 , and in the range [145 − 160] µm for φ 0 = 3.57 × 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 φ 0 . Equation 1 has been derived by a scaling theory that involves δ 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 δ 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 δ 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 The main caveats for this formulation are that it strongly depends on having a correct scaling for δ 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.

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 (t) and consequently for the mass of particles accumulated at the bottom of the tank m b (t). 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 (t) and m b (t) 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 (t) 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 + t a where t a is the time during which the PBL initially grows beneath the interface at the individual particle settling velocity, i.e., t a = δ PBL /V s , and t a 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 = (h 2 − δ PBL ) /V f . We therefore see that t a strongly depends on δ 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.

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;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; 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 ∼ 10 4 in our configuration. This value differs from the one suggested by analogy with thermal convection (Gr c ∼ 10 3 ) . 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 δ 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 200021_169463.

Supplementary Material
June 16, 2021 1. DESCRIPTION OF THE FINITE DIFFERENCE SCHEMES 1.1. First-order upwind finite difference scheme In the following description, we describe only the one-dimensional scheme as the method is easily generalised to higher dimensions by simply applying the procedure separately to each dimension (Ferziger and Peric, 2002). The finite difference method is based on the approximation of the derivatives at the node locations of a discretised domain and can be applied to both uniform and non-uniform meshes. However, we describe here the case where the numerical domain is uniformly discretised with the spatial step δ x, such that the domain is divided in a set of equally-spaced points {x 0 , x 1 , ..., x i ..., x n }. We also consider the one-dimensional conservation equation where a is the transported information and u the advection velocity. In our discrete domain, we use the notation a n i = a (x = x i ,t = t n ),with n denoting the time coordinate and i the spatial coordinate. Among the numerous methods used to approximate the derivative ∂ a/∂ x at location x i , the Taylor expansion is the most common. Following this procedure, there are two ways to estimate the derivative with a first order accuracy: • Using the Taylor expansion at the location x i+1 = x i + δ x, we get the forward difference ∂ a ∂ x F ≈ a n i+1 − a n i δ x .
• Using the Taylor expansion at the location The central finite difference approximation can be determined by combining the forward and backward differences. Thus, we have the second order accurate discretisation The upwind finite difference method scheme is an adaptive procedure to discretise the problem based on the direction of propagation of the information. The estimation of the quantity a n+1 i depends on the sign of u: • If u > 0, the backward difference is used a n+1 i = a n i − δt δ x u a n i − a n i−1 .
Modelling gravitational instabilities • If u < 0, the forward difference is used a n+1 i = a n i − δt δ x u a n i+1 − a n i .
All these properties of the first order upwind scheme guarantee a strong stability of the solution providing the Courant-Friedrichs-Lewy (CFL) condition is satisfied (Courant et al., 1928) δt δ x u ≤ 1.
1.2. Third-order Weighted Essentially Non-Oscillatory (WENO) finite-difference scheme In the following description, the numerical domain is discretized using the same set of points as presented in the previous section. The low-diffusive WENO scheme belongs to a family of high-resolution methods and was developed in order to solve hyperbolic partial differential equations of the form where, in our case, the flux takes the form h (a) = ua (x) . The WENO scheme provides a third order accurate method in smooth regions, i.e., where the spatial gradient is small. On the other hand, the scheme adapts where the gradient is high. In these cases, the accuracy tends toward second-order. This adaptive aspect of the WENO procedure ensures the suppression of spurious oscillations, predicted by the Godunov theorem , around shocks (regions of high gradient). The principle is to build a convex combination of interpolants for the flux at given points of the domain, using different stencils. The third-order WENO procedure uses two adjacent stencils of two points each ( Figure  SS1). Figure S1: Numerical stencil used for the 3rd order WENO procedure It is possible to approximate the spatial derivative with the commonly used half-node flux: where h i+ 1 2 = h a x i+ 1 2 is the numerical flux at the half-node location x i+ 1 2 (i.e., the central position between the points x i and x i+1 ). There are three ways to construct the half node flux h i+ 1 2 . Firstly, we can use polynomial interpolants of degree 1 p 1 (x) and p 2 (x), respectively, in the two-point stencils 1 and 2. Then, the flux at the location x i+ 1 2 is given either by h and Note that the flux approximations for each stencil in this case have an accuracy of second order. The last way to determine the flux h i+ 1 2 is with an interpolating polynomial p 3 (x) of degree 2 inside the three-point stencil 3 (i.e., the union of stencils 1 and 2 in Figure SS1). Then, the third-order flux at the location x i+ 1 2 is given by The three interpolations of the half-node flux presented above are efficient assuming that the function h is smooth through the associated stencil and it is even possible to write the third-order flux h where γ 1 = 3/4 and γ 2 = 1/4. However, the presence of any discontinuity would break the stability of the procedure, introducing spurious oscillations in the calculated solution. The treatment of this aspect constitutes the essence of the WENO method which retains the property of relating the total flux h i+ 1 2 to a convex combination h (1) i+ 1 2 and h (2) i+ 1 2 , similarly to equation S13, but including non-linear weights where ω 1 and ω 2 are functions of the smoothness of h and must satisfy ω 1 + ω 2 = 1. If h is smooth in all the stencils, ω i → γ i and ensures a third-order accuracy. Conversely, the presence of any discontinuity in the i − th stencil means ω i → 0, decreasing the accuracy to second-order. This property is guaranteed by determining ω i and, in the third-order WENO case, α 1 = 2/3, α 2 = 1/3, β 1 = (h i+1 − h i ) 2 and β 2 = (h i − h i−1 ) 2 (β j are referred to as smoothness indicators). A positive mathematical coefficient ε is introduced in order to avoid any division by zero for the calculation ofω j . The WENO reconstruction presented above is valid for u > 0. In order to include potential changes of the flow direction, we use an upwind splitting of the flux. Thus, as a mirror image of the procedure described above, we have for u < 0 and With an appropriate time discretisation, the third-order WENO procedure remains a Total Variation Diminishing (TVD) scheme, i.e., ∑ i |a n+1 i+1 − a n+1 i | ≤ ∑ i |a n i+1 − a n i | (Harten, 1983). This property avoids the introduction of new local extrema in the solution. Therefore, a TVD third-order WENO scheme is given using an iterative third-order Runge-Kutta method for the time discretisation: a (1) = a n + δtW (a n ) , and a n+1 = 1 3 a n + 2 3 with W (a) is the spatial derivative given by the WENO algorithm described above.

EXTENSION OF ANALYTICAL MODEL OF HOYAL ET AL. (1999)
2.1. Particle concentration For a particle suspension placed above a denser fluid, an estimate for the evolution of the particle concentration in the lower layer has been derived by . We provide here an extension of this formulation.

Upper layer
First, let's consider a quiescent upper layer of thickness h 1 with an initial particle concentration C 1 (0). For a quiescent settling at the velocity V s , the constant flux of particles leaving the upper layer is given by −AV s C 1 (0). Thus, the evolution for the mass of particles is described by the equation: Assuming that the mass of particles depends on the concentration as: A being the horizontal cross section of the tank (where i = 1 for the upper layer and i = 2 for the lower layer), the temporal evolution of the particle concentration in the upper layer is the solution of: Then:

Lower layer
Now, let's defineṀ in as the flux of particle entering the lower layer (i.e. the particles arriving from the upper layer) andṀ out the flux of particles leaving (i.e. particles that deposit at the bottom of the tank). Thus, the mass of particle in the lower layer M 2 (t) is governed by the equation • There is a time t a before which particles have not yet reached the bottom of the tank and evidentlẏ M out = 0. Still assuming a quiescent upper layer, we have also: So, combining S25, S28 and S29 the particle concentration in the lower layer for t < t a is the solution of : that is to say (assuming C 2 (0) = 0): • Once the first particles have reached the bottom of the tank (i.e. t ≥ t a ), the flux of particles leaving the lower layer becomesṀ out = AV s C 2 (t) .
Similarly to the previous point, combining equations S25 (which assumes that the lower layer is turbulently convecting and is homogeneous), S28, S29 and S31 we obtain Using the initial condition C 2 (t a ) = (V s /h 2 )C 1 (0)t a (continuity with S31), the solution of this equation is: It is evident that when t a = 0, S34 is equivalent to the original formulation of  (equation 4 in the Introduction). • However, S34 is only valid for the condition that particles keep settling across the interface i.e. for a time t < t lim (where t lim = h 1 /V s ). We have extended the formulation to later times once all particles from the upper layer have settled across the interface i.e. for a time t ≥ t lim . As there are no longer particles in the upper layer, the fluxṀ in drops to zero. Two possibilities are now available. On one hand, if we consider that the lower layer is still turbulently convecting after t lim , equation S33 becomes: Assuming that all particles have crossed the interface by the time t lim , we substitute t lim into S34 in order to find the initial condition for equation S35. Thus the initial condition is now: C 2 (t lim ) = C 1 (0) 1 + V s h 2 t a − 1 e A Area of the plane through which the particle mass flux is computed b Coefficient used to compute the force term in the  formulation c i Local particle velocity associated to a given lattice type in the LBM c s Speed of sound used in the LBM C 1 Particle concentration in the upper layer C 2 Particle concentration in the lower layer D, D c , D S Diffusion coefficients respectively for the density-altering quantity, the particle field and the sugar D p Particle diameter D z Differential operator e z Vertical unit vector f i , f eq i Particle population and equilibrium distribution function F Body force term Mass of particles in the lower layeṙ M in ,Ṁ out Particle mass flux entering and leaving the lower layer respectively NE Numerical diffusion error in the first order finite-difference scheme p, p c Fluid pressure and characteristic pressure in the viscous scaling q Power used in the law defining V f as a function of φ S, S 0 , S * ,S,Ŝ Sugar concentration, sugar concentration at time t = 0 and S * (= S 0 ) is used in the nondimensionalisation, base state and perturbation amplitude for the sugar concentration Sc i Schmidt number St Stokes number t,t a ,t c ,t lim Time, time for the particle first arrival at the bottom of the tank, characteristic time in the viscous scaling and time at which all particle have settle across the initial interface. T Time at which the instability starts growing u, u f Transport velocity and fluid velocity U Characteristic velocity V f ,V s Respectively the fingers velocity and individual particle settling velocity W Right term in the matrix form of the eigenvalue problem in the LSA x = (x, y, z) Position vector and associated 3D components z φ , z s Error function parameters used for the base states of volmue fraction and sugar respectively α Sugar expansion coefficient Γ H , Γ H i Power spectral density (PSD). The i index stands for the initial PSD δt, δ x Temporal and spatial steps δ PBL PBL thickness ∆t Integration time for the particle mass flux η Power used in the law defining V f as a function of D p µ Fluid dynamic viscosity ν Fluid kinematic viscosity ρ, ρ 0 , ρ f , ρ p , ρ PBL , ρ blk Density, fresh water density, fluid density (including sugar), particle density, PBL density, bulk density (including sugar and particles) σ , σ sim , σ LSA,i Growth rate of the instability, growth rate computed from simulations and growth rate predicted by LSA τ Relaxation coefficient in the LBM ϕ,φ, ϕ ,φ Arbitrary variable, associated base state, perturbation and perturbation amplitude φ , φ tot , φ i , φ 0 , φ * ,φ ,φ Particle volume fraction, total particle volume fraction (polydisperse case), volume fraction of the i − th size class, initial particle volume fraction, φ * (= φ 0 ) is used in the nondimensionalisation, volume fraction base state, perturbation amplitude ψ,ψ,ψ Stream function, associated base state, perturbation amplitude ω,ω,ω Vorticity, associated base state, perturbation amplitude s Table S1: List of symbols used in the main manuscript