The experiments are performed in an experimental water tank 30.1 × 30.1 × 50.0 cm³ that is divided into two layers by a horizontal barrier (Manzella et al., 2015; Scollo et al., 2017; Figure 2; H 1 = 5–21 cm and H 2 = 25 cm are the thicknesses of the upper and the lower layers, respectively). Initially, the density profile is stable. The upper layer (ash cloud analogue) is a mixture of particles and fresh water, with initial particle concentrations C u ranging from 0.0007 to 10 g l⁻¹ (mixture density from 997 to 1,007 kg m⁻³), and the lower layer (atmosphere analogue) is a denser sugar solution kept at a constant density of 1,008 kg m⁻³ (sugar concentration of 29 g l⁻¹). Experiments are performed at ambient temperature. This experimental configuration is inspired by the classical work of Hoyal et al. (1999b) on settling-driven gravitational instabilities and experiments on Rayleigh-Taylor instabilities (Linden and Redondo, 1991; Dalziel, 1993). This minimalist design allows us to study the effect of particle concentration on settling-driven gravitational instabilities and finger dynamics in isolation from other processes and variables.

The particles are spherical glass beads, which are suitable for modeling the behaviour of natural ash (Manzella et al., 2015; Scollo et al., 2017). They have a median diameter D 50 of 41.5 µm and a sorting σ = ( D 84 − D 16 ) / 2 (Inman, 1952) of 14 µm (see Supplementary Figure S1), as measured by laser diffraction using a BetterSizer S3 Plus. The particle density is measured, using Helium pycnometery with a Ultrapyc 1200e, to be 2,519.2 ± 0.1 kg m⁻³.

Particles are kept in suspension by continuously mixing the upper layer top to bottom for 20–25 s with an agitator composed of a 23 × 7 cm² millimetric mesh, that is large enough to mix the entire upper layer and produce homogeneous and repeatable mixing. Mixing is stopped 5 s before the experiment begins. The separation between the two layers is then removed manually by sliding the barrier out of its slot in 0.9–1.3 s. A major challenge of this experimental configuration is that mixing is generated across the density interface when removing the barrier and is clearly visible in the first 10 s of the experiments, affecting the early development of settling-driven gravitational instabilities. The perturbation is particularly strong near the back wall of the tank, so we reduced the length of the tank to 7.5 cm with rigid vertical separators to attenuate the effect of the vorticity (Manzella et al., 2015; Scollo et al., 2017). After barrier removal, particles start settling into the lower partition of the tank (see Supplementary Video S1). A continuous Nd:YAG planar laser (Genesis CX-SLM by Coherent) with a wavelength of 532 nm illuminates the experiments from the side of the tank and a sCMOS camera (HiSense Zyla by Dantec Dynamics), with 16 bit color depth, is used to capture images of the experiments at 10 Hz with image dimensions of 1,500 × 2,100 pixels, resulting in a resolution of 55 px mm⁻¹. We find that the sugar concentration is sufficiently small that it has negligible effect on the light intensity received by the camera.

Gravitational instabilities in our setup can theoretically develop through two different mechanisms, double diffusion and particle settling (Hoyal et al., 1999a; Hoyal et al., 1999b; Burns and Meiburg, 2012). Double diffusive effects arise when two density-altering fluid properties (i.e., sugar and particles in the lower and the upper layers, respectively) diffuse at differential rates. The faster diffusion of one component relative to the other can lead to a local increase of the bulk density at the interface that triggers gravitational instabilities. The rate at which the bulk density increases in double-diffusive systems is therefore controlled by the diffusion of the fastest diffusing component. Conversely, settling-driven gravitational instabilities are generated by particles settling across the interface, causing the upper part of the lower layer (i.e., the PBL) to become heavier than the fluid below. Gravitational instabilities produced by particle settling are therefore controlled by both the vertical velocity of the particles and the initial particle concentration of the upper layer. We determine the mechanism that provokes the formation of instabilities in our experiments by calculating the ratio of the double diffusive and settling fluxes, F D and F I , respectively, at the density interface (Green, 1987; Hoyal et al., 1999a; Carazzo and Jellinek, 2013) F ∗ = F D F I . (2)

Replacing the fluxes by their expressions F D = 1 20 ρ u ( g κ f β m ) 1 3 ( C u / ρ p ) 4 3 and F I = C u g d 2 ( ρ p − ρ f ) / 18 μ (Hoyal et al., 1999a; Carazzo and Jellinek, 2013), where κ f is the diffusion coefficient of the fastest diffusing substance, μ the dynamic viscosity of the fluid, C u the initial particle concentration in the upper layer (in mass per unit volume), ρ p the particle density, d the particle diameter, β m ∼ ( ρ p − ρ f ) / ρ p the volumetric expansion coefficient of a particle suspension and ρ u and ρ f the density of the particle suspension and of the fluid phase, respectively, Eq. 2 can be written as F ∗ ≈ 9 μ ρ u ( κ f β m C u ) 1 3 10 d 2 ( ρ p − ρ f ) ( ρ p 2 g ) 2 3 . (3) F ∗ >> 1 means that double diffusion dominates the mass flux across the interface, whilst particle settling dominates for F ∗ << 1. In our experiment, the particles are the fastest diffusing substance, with a maximum hydrodynamic diffusion coefficient κ p of 1.2 × 10⁻⁷ m² s⁻¹ (Lee et al., 1992; Martin et al., 1994). We find that F ∗ is between 8 × 10⁻⁴ and 2 × 10⁻² in our experiments, suggesting that settling is the principal mechanism provoking the formation of fingers, with a negligible contribution of double diffusive effects.

Our experiments are dedicated to correctly reproducing the processes affecting ash sedimentation beneath volcanic clouds in a small, simplified configuration. The difference in complexity and scale between the natural phenomenon and small-scale experiments raises the problem of the applicability of the analogue experimental results. To address this, we perform a scaling analysis (Burgisser and Bergantz, 2002; Burgisser et al., 2005; Carazzo and Jellinek, 2012; Kavanagh et al., 2018; Roche and Carazzo, 2019). Dimensionless numbers relevant to our problem are calculated for volcanic ash clouds and compared with particle suspensions (Figure 2D). We also compare the dynamical regimes of natural ash and experimental fingers based on dimensionless numbers (Figure 2E). The values of the parameters associated with volcanic clouds and fingers are obtained from the literature (Table 1; Burgisser et al., 2005; Carazzo and Jellinek, 2012; Manzella et al., 2015).

The Reynolds number R e = ρ V L / μ , with μ the fluid dynamic viscosity and V and L the characteristic velocity and length scale of the flow, respectively, characterises the flow behavior by comparing the inertial to viscous forces in the fluid. R e is greater in natural volcanic clouds and fingers, where it is far above the mixing transition ( R e > 1,000–4,000) and is fully turbulent. Hence, we can expect particle collisions and entrainment of the ambient fluid to be enhanced in natural ash fingers, compared with their experimental counterparts. However, the flows in both experiments and volcanic clouds are inertia-controlled, suggesting that they are comparable despite the fact that velocity fluctuations (i.e., turbulence) are greater in natural flows.

In situations involving particle settling in fluids, the particle Reynolds number R e p = ρ f V p d / μ is used to assess the properties of the flow surrounding particles and the subsequent drag force acting on them. d is the particle diameter and V p is the individual particle settling velocity given by the Stokes terminal velocity in the experiments V p = g d 2 ( ρ p − ρ f ) 18 ρ f ν . (4) At subcritical R e p < 3 × 10 5 , the drag factor f , which assesses the importance of the fluid resistance exerted on the particles, can be expressed as a function of R e p and two drag correction coefficients depending on the particle shape k N and k S by Bagheri and Bonadonna (2016a), Bagheri and Bonadonna (2016b) f = 24 k S R e p [ 1 + 0.125 ( R e p k N / k S ) 2 3 ] + 0.46 k N 1 + 5330 R e p k N / k S . (5) For multiphase flows involving particles, the Stokes S t and Sedimentation Σ numbers quantify the momentum transfer between the fluid phase and the particles (Burgisser et al., 2005; Roche and Carazzo, 2019). They are calculated as S t = ( ρ p − ρ f ) d 2 18 µ f ( Δ U δ f ) ( 1 + ρ f 2 ρ p ) , (6) and Σ = ( ρ p − ρ f ) d 2 18 µ f ( g Δ U ) = V p Δ U , (7) with Δ U the characteristic velocity fluctuation over a characteristic distance δ f . These two parameters allow us to assess the coupling between particles and fluid. When S t << 1 and Σ << 1, particles are strongly coupled with the fluid. This is the case for both natural ash and experimental fingers, satisfying the assumption that particles remain coupled with the fluid for fingers to form (Hoyal et al., 1999b; Carazzo and Jellinek, 2012).

Previous studies have invoked a convective mechanism for finger formation whose initiation depends on the Grashof number G r , which can be estimated using a combination of experimental and field observations (see section Finger Length Scales in Experiments: PBL Thickness, Finger Width and Finger Spacing). Our analysis shows that G r is much greater in the natural systems than in the experiments but that, in both configurations, it exceeds the critical Grashof number G r c = 10 3 (Hoyal et al., 1999b) for the development of settling-driven gravitational instabilities.

For scaling the density ratio between the particle-laden layer and the underlying fluid, we calculate the Atwood number A t = ( ρ − ρ a ) / ( ρ + ρ a ) . A t is commonly used to parameterise the density difference between two fluid layers, notably in experiments on Rayleigh-Taylor instabilities (Dalziel, 1993; Wilson and Andrews, 2002). Here, A t << 1 in all configurations, showing that the difference in density between the two layers is very small in nature as well as in experiments and that the instability develops from small initial density differences in both situations.

For settling-driven gravitational instabilities, the natural and experimental ranges of R e p , A t , S t and Σ all overlap while they are different for R e and G r , but above critical threshold values that ensure that experiments are comparable to the volcanic phenomenon. We also note that the dimensionless numbers systematically cover a wider range in nature than in experiments (Table 1). This can firstly be explained by the wide variety of eruptive source and atmospheric conditions which means eruptive clouds can be associated with a wide range of particle diameters, characteristic velocities and length scales. Moreover, a second explanation is that descriptions of ash fingers remain rare, meaning they are currently poorly constrained, with high associated uncertainties. Altogether, this scaling analysis shows that our experiments do a good job reproducing most dimensionless numbers associated with natural clouds and fingers, although the variability of natural phenomena means that they are associated with larger ranges of R e and G r , potentially extending to greater values than in our experiments. However, it is inevitable that laboratory models of volcanic clouds cannot capture the full range of variability of the natural system (Carazzo and Jellinek, 2012; Kavanagh et al., 2018; Roche and Carazzo, 2019). Despite this, given that that the ranges of R e and G r in both the experiments and natural clouds are close to or above expected transitional values (Reynolds, 1883; Hoyal et al., 1999b), we can regard our experiments as suitable analogues to study settling-driven gravitational instabilities at the base of volcanic clouds.

In order to independently obtain both the particle concentration C and fluid density fields, we imaged the particles using a set-up derived from Particle Image Velocimetry (PIV; Keane and Adrian, 1992; Grant, 1997; Adrian, 2005), and the fluid phase using Planar Laser Induced Fluorescence (PLIF; Koochesfahani, 1984; Crimaldi, 2008), in separate repeated experiments. These two techniques have previously been applied separately or simultaneously (Borg et al., 2001; Hu et al., 2002) to assess the particle concentration and velocity fields, as well as fluid properties such as temperature and density, in other aqueous experiments. Here we apply particle imaging and PLIF on separate experiments where we kept the experimental set-up in the same configuration, including identical starting conditions, i.e., particle size and concentration and fluid density contrast, as well as identical imaging conditions, i.e., camera and laser positions and settings. Ensuring the same experimental and imaging conditions is critical in order to enable the combination of the results given by the two techniques (Borg et al., 2001). Moreover, we ensured that results from PLIF experiments were reproducible before combining with results from particle-imaging experiments (Supplementary Figure S2).

We developed a calibration that links the particle concentration C to the light intensity scattered by particles and received by the camera sensor, with the assumption that the light intensity observed in the experiments is linearly related to the particle concentration. During calibrations, the water tank is first filled with a fluid of uniform and known particle concentration (ranging from 0 to 6 g l⁻¹ at 0.5 g l⁻¹ concentration steps). The tank is then illuminated with different laser powers of 0.50, 0.45, 0.40, and 0.35 W for each particle concentration. In order to erase short temporal variations in particle concentration, calibration images are acquired at 0.1 s intervals and then averaged over 10 s. Finally, we perform an individual pixel-by-pixel calibration relating pixel intensity to particle concentration for each laser power. The pixel digital level (i.e., received light intensity) is positively linearly correlated with the particle concentration and the coefficient of determination is high, R ² = 0.95 (Figure 3A). However, the light intensity through the width of the experimental tank is dependent on the particle concentration and exponentially attenuates with x, the distance from the tank wall closest to the laser. Hence, the quality of the linear regression between light intensity and particle concentration diminishes with distance from the light source (Figures 3B,C). For accurate particle concentration measurements, we therefore calculate concentration only in pixels with a coefficient of determination R ² > 0.95. This corresponds to the pixels located in the 12 cm of the tank closest to the laser, limiting the particle concentration measurements to approximately one third of the tank width. We assume that finger characteristics in this section are representative of those throughout the whole tank. Particle concentration measurements have already been employed to characterize settling-driven gravitational instabilities by Hoyal et al. (1999b) and Manzella et al. (2015). In both studies, the authors placed light sources behind or beside their experimental tanks and related light attenuation to the depth-averaged particle concentration in interrogation areas. With the present particle imaging scheme, we managed to improve the resolution of the concentration measurements down to the pixel-scale, which corresponds to a resolution of 0.18 mm in the horizontal and vertical directions. Furthermore, our measurements within the laser plane differ from depth-averaged particle concentration measurements and offer a visualization of the 2D particle concentration. We compare our particle concentration measurements obtained using this novel calibration procedure against those of Manzella et al. (2015). The two measurement techniques are in good qualitative agreement (Supplementary Figure S3). It is worth noting that whilst particle concentration measurements in the lower layer, into which the laser light enters directly, are highly accurate, corresponding measurements in the upper layer are affected by the refractive index change between the sugar solution and the particle suspension.

Understanding the behavior of the fluid phase associated with the occurrence of settling-driven gravitational instabilities is crucial to improving our general comprehension of the mechanisms driving the instability. In particular, it is important to quantify the contribution of the fluid phase on the generation of an unstable density profile associated with settling-driven gravitational instabilities. For this purpose, we applied the PLIF technique to a series of separate experiments. The PLIF approach yields quantitative, non-intrusive, measurements of fluid concentrations at high resolution (Koochesfahani, 1984; Crimaldi, 2008). It involves using a dye that generates fluorescence with an intensity linearly proportional to its concentration D when illuminated by a laser. In experiments on fluid mixing, by adding dye to one endmember, the fluorescence intensity can be used as a tracer, with the fluorescent dye concentration related to the proportion of dyed fluid (Linden and Redondo, 1991; Troy and Koseff, 2005; Dossmann et al., 2016). In the case presented here, mixing occurs between fluids of different density. The bulk density of the fluid phase thus changes during mixing and is quantified by the fluorescent dye concentration.

During our PLIF experiments, the upper layer is doped with Rhodamine 6G (R6G), a fluorescent dye with an absorption peak at 530 nm (Zehentbauer et al., 2014), close to the laser wavelength (532 nm), for a maximum absorption efficiency. The emission spectrum of R6G ranges from 510 to 710 nm with a peak at 560–580 nm, depending on the solvent (Zehentbauer et al., 2014). We isolate the R6G fluorescence from other light sources with a smaller wavelength (i.e., primary laser emission, particle scattering) by equipping the camera with a high-pass filter at 570 nm. This procedure allows us to image the R6G distribution only, as a proxy for the concentration of the upper fluid phase.

We perform a calibration in order to relate the local R6G dye concentration D to the digital level of individual pixels (i.e., fluorescence intensity). In PLIF measurements, the digital level F increases linearly with the R6G dye concentration D (Borg et al., 2001; Crimaldi, 2008) F = ( ϕ ϵ P a ( r , θ ) γ ( r , θ ) Δ A ) D + R , (8) in which Δ A is the pixel area, a ( r , θ ) the attenuation along the path, ϵ the dye absorption coefficient, P the laser power, ϕ the dye quantum efficiency,   γ ( r , θ ) the spatial intensity distribution, R the residual light not related to R6G fluorescence, and r and θ the radial and the angular components, respectively, of the 2D polar coordinates describing the position of points within the laser plane. Keeping the same particle concentration, fluorescent dye, laser properties, camera settings and experimental configuration in all calibrations and experiments, Δ A , a ( r , θ ) , P , ϕ , ϵ , γ ( r , θ ) and R can reasonably be assumed to be constant. Equation 8 thus reduces to F = p D + R , (9) where p and R are constants that can be determined by linear regression between F and D .

However, changing the particle concentration modifies the spatial intensity distribution because the particles have an attenuation effect on a ( r , θ ) and, therefore, the determination of p . In fact, complex optical effects including light scattering and shadowing take place in the presence of particles and the light intensity I decreases following approximately exponential curves, whose slope steepens with particle concentration (Figure 4A). Hence, in order to study the mixing between a particle-laden and a clear fluid, it is necessary to account for this effect (Borg et al., 2001). Whilst calibrations are performed with a uniform particle concentration, the spatial distribution of particles in the experiments can be very complex. Therefore, since the particle and R6G concentrations cannot be assessed simultaneously in this experimental set-up, we conduct PLIF calibrations using different particle concentrations. The difference between calibrations allows us to assess the uncertainty on measurements where the particle concentration is unknown. To obtain minimum uncertainties, we perform calibrations for particle concentrations of 0, 1, and 2 g l⁻¹, where the light intensity diminishes to a minimum of 33% of its initial value (Figure 4A). The effect of particle concentration on R6G fluorescence through the width of the tank is shown in Figure 4B. It is clear that small particle concentration differences can drastically reduce the intensity of fluorescence by attenuating the spatial light distribution inside the tank.

We conduct PLIF calibrations varying the particle and R6G concentration whilst keeping a constant laser power of 1.15 W. The water tank is filled with a uniform dye concentration ranging from 0 to 8 μg l⁻¹ for three particle concentrations: 0, 1, and 2 g l⁻¹. Calibration images are averaged over 10 s in order to erase short temporal variations in fluorescence intensity. For each pixel, we fit the digital level and the R6G concentration to Eq. 8 to determine p and R . In contrast with the particle concentration calibrations, the quality of the linear regression is good throughout the tank for a uniform particle concentration, where a ( r , θ ) is constant. However, the slope p is affected by the presence of particles as the fluorescence intensity diminishes with increasing particle concentration (Figure 4C), yielding an uncertainty on the measurements of R6G concentration and ultimately on the fluid density measurements. In the region of particular interest where the PBL forms, located within the first 5 cm below the barrier, we estimate the average uncertainty of the final fluid density measurement to be 0.1 and 0.8 kg m⁻³ for 1 and 2 g l⁻¹ PLIF experiments, respectively (Supplementary Figure S4).

We performed two types of experiments, with different objectives. Type A experiments are each repeated at least three times with the aim of separately characterising the fingers at different particle concentrations, and the density of the fluid phase (Table 2). In the first subset of type A experiments (A1–A10), we image the particles and vary the particle concentration from 1 to 10 g l⁻¹. In the second subset of type A experiments (A11–12), we repeat the same experiments but instead measure fluid phase properties using the PLIF technique, only with particle concentrations of 1 and 2 g l⁻¹. Type A experiments have a constant upper layer thickness of 13.5 cm and are imaged for 90 to 120 s after barrier removal.

In type B experiments, we explore the conditions that favor the formation of settling-driven gravitational instabilities by varying both the upper layer thickness (from 5 to 21 cm) and the particle concentration (from 7 × 10⁻⁴ to 3 g l⁻¹; Table 2). The maximum upper layer thickness is limited to 21 cm by the size of the water tank and we select the minimum thickness to be 5 cm in order to keep the upper layer much greater than the size of perturbations that could be induced by barrier removal. All type B experiments are particle-imaging experiments used to identify the presence or absence of settling-driven gravitational instabilities and are imaged for at least 180 s after barrier removal.

For PLIF experiments (A11–12), we calculate the local concentration of upper fluid X f that corresponds to the proportion of upper layer fluid (containing the fluorescent dye) within the fluid phase as X f = D / D u , with D the local dye concentration and D u the initial R6G concentration of the upper layer. We then determine the fluid density ρ f from the concentration of upper fluid as ρ f = X f ρ w + ( 1 − X f ) ρ s , (10) where ρ w is the density of fresh water (that initially forms the fluid phase of the upper layer) and ρ s is the density of the sugar solution in the lower layer. Moreover, PLIF experiments provide reproducible results, with fluid density profiles very similar from one experiment to another. This allows us to reasonably combine the fluid phase density with the particle concentration, obtained in experiments A1–10, in order to calculate the bulk mixture density ρ m as ρ m = X p ρ p + ( 1 − X p ) ρ f , (11) with ρ p the particle density and X p = C / ρ p the particle volume fraction calculated from the measured particle concentration C . At the density interface, the reduced gravity is estimated as a function of ρ m and ρ s as g ′ = g ρ m − ρ s ρ s . (12) For the PBL of density ρ , the reduced gravity is calculated as g ′ = g ( ρ − ρ s ) / ρ s .

Whilst PLIF experiments used to infer the fluid phase density involve only small particle concentrations of 1 and 2 g l⁻¹, we additionally calculate the bulk density in experiments involving higher particle concentrations by estimating the fluid phase density from C u = 2 g l⁻¹ PLIF experiments. To do so, we assume that the spatial distribution of the fluid phase is only weakly affected by the particle concentration. Given that the particle concentration can reach up to 10 g l⁻¹, this is a strong assumption that we can only test for C u = 1 and 2 g l⁻¹ PLIF experiments for which fluid phase density profiles are indeed very similar (Supplementary Figure S2).

Type A particle-imaging experiments (A1–10) are processed in order to quantify the effect of particle concentration on the dynamic and geometric properties of fingers. The finger velocity is obtained by manually tracking the front position (i.e., lowest position of the interface between the finger and the lower layer fluid) over time. We measure the PBL thickness throughout experiments using both visual observations and particle concentration profiles. At the lower boundary of the PBL there is an abrupt change in particle concentration that corresponds to the transition from the particle suspension to the particle-free lower layer. We therefore define the PBL base as the position where the vertical gradient of the particle concentration reaches its minimum value. Results show that the position of the lower PBL boundary (and therefore the PBL thickness) reaches a constant height after 10 to 20 s (Supplementary Figure S5), when the initial effect of barrier removal disappears. The PBL thickness is finally defined as the vertical separation between the initial density interface, represented by the barrier, and the lower PBL boundary (Figure 5). When fingers develop, they form distinct particle-rich columns that descend into the lower layer and are characterised by bright regions in the experiments. Along a horizontal transect crossed by fingers, the light intensity profile has a peak at each finger location (Supplementary Figure S6). We therefore determine the number of fingers by counting the most prominent peaks along a transect located 4 cm under the upper layer, where fingers are clearly developed. This counting procedure is sensitive to detection parameters such as the vertical position of the transect or threshold values for peak prominence, width, and separation. Therefore, in order to obtain repeatable results comparable to other studies, these parameters have been selected to reproduce manual finger detection (Scollo et al., 2017) within ±1 finger. In contrast to previous methods, this allows us to quantify the temporal evolution of the number of fingers throughout the experiment duration. We assess the finger spacing λ and width W from visual observations (Figure 5). These geometrical parameters are measured 10 to 20 s after the experimental onset and at the same time as the PBL thickness, when the effect of barrier removal becomes small. λ corresponds to the distance separating two fingers and is measured at the base of the PBL, whilst W is measured at the head of the fingers (i.e., their thickest point). The uncertainty is quantified by the standard deviation from measurements of multiple fingers in each experiment. Key experimental measurements are summarised in Table 3.

Type B experiments are exclusively used to determine the presence or absence of fingers over a wider range of initial conditions than type A experiments. The presence of fingers is determined by visual identification of i) downward-moving particle-laden plumes, with columnar or bulbous shapes, that form at the initial density interface and inside which particles move faster than individually and ii) the tendency of fingers to initiate and maintain layer-scale convection in the lower partition of the tank (Hoyal et al., 1999b). We do not measure the evolution of particle concentration in type B experiments, for which the calibration presented in Particle Imaging. does not apply because of the low initial particle concentrations employed.

Particles are initially contained in the upper part of the tank and they start sedimenting through the initial density interface after barrier removal. Subsequently, the average particle concentration decreases over time in the upper layer whilst it increases in the lower layer. In all experiments, we observe that the particle concentration in the PBL does not exceed that of the upper layer. Horizontally averaging the particle concentration at each time, we obtain particle concentration profiles for experiments A1–10. Particle concentration profiles are then combined with fluid phase density profiles (Eq. 10) to approximate the bulk mixture density profile (Eq. 11).

Figure 6A shows the evolution of bulk mixture density with time for experiment A8 (8 g l⁻¹; Table 2), chosen as an example since the relatively high particle concentration leads to a pronounced unstable density profile. Just after barrier removal ( t = 3 s), the bulk density increases with depth in the experimental apparatus and the density profile is stable. However, particles sink below the interface and mix with the dense lower layer fluid. Consequently, an unstable heavy region associated with the PBL grows in thickness and density below the upper layer because of the added particle concentration. For experiment A8, this layer is particularly pronounced after the time t = 12 s, where it exceeds the density of the lower layer by about 1.3 kg m⁻³.

An initial increase in density in the region immediately below the initial interface is observed in all experiments. Figure 6B shows the maximum excess density (i.e., difference between the density of the unstable region and that of the lower layer) vs. initial upper layer particle concentration, for experiments A1–A10 (Table 2). The maximum excess density gradually increases from 0.15 kg m⁻³ for a particle concentration of 1 g l⁻¹ to 1.4 g m⁻³ for a particle concentration of 10 g l⁻¹.

The perturbation to the system caused by the barrier removal at the experiment onset prevents us identifying the timescale of growth of the PBL thickness, which is expected to be δ / V p (Martin and Nokes, 1989). Because of this experimental limitation, we therefore quantify the timescale at which the density increases in the PBL by measuring t i , the time to develop an unstable layer with excess density ≥0.15 kg m⁻³. We find that t i decreases with C (Supplementary Figure S7; Supplementary Videos S2, S3), with t i = 40–60 s at C u = 1 g l⁻¹ and t i = 4–6 s at C u = 10 g l⁻¹ (Figure 6C). Measurements are well fitted by the relation t i = ( α / g ′ ) , where α = 0.2 m s⁻¹ is an empirically-fitted constant, dimensionally homogeneous to velocity. By analogy with the expected timescale of PBL growth, we presume that the constant α , which governs the increase of the density inside the PBL, is inversely proportional to the individual settling speed of particles.

Unlike particles, the upper layer fluid phase is not significantly entrained inside the lower layer after barrier removal and fluid mixing between the upper and the lower layer is limited in proportion and vertical extent (Figure 7). However, we notice the development of a thin 2–3 cm mixing zone below the interface with intermediate concentrations of upper fluid (10–60%). This zone is affected by disturbances created by the initial barrier removal. Below this zone, the spatial distribution of the concentration of upper fluid shows that only small amounts of upper fluid are entrained within fingers. This entrained fluid appears as thin leaks of upper fluid [similar to those observed by Parsons et al. (2001)] with concentrations of upper fluid up to 6–7% (Figure 7B). Once entrained downward, portions of the upper layer are seen to buoyantly rise back up once fingers reach the bottom of the tank (Supplementary Videos S4, S5). On average, the concentration of upper fluid in the lower layer increases with time but remains low (<3%) (Figure 7C). There is an initially rapid increase at the start of the experiment due to the removal of the barrier. This is then followed by small oscillations that are linked with the thickening and retraction of the PBL at the location where fingers are formed. This thickening induces a perturbation of the lower density interface that results in an increase of the average concentration of upper fluid. The PBL later retracts to its initial position after detachment of the fingers, resulting in a decrease in the average concentration of upper fluid.

The average finger velocity increases with particle concentration (Figure 8), from 0.35 cm s⁻¹ at C u = 1 g l⁻¹ to 1.4 cm s⁻¹ for C u = 10 g l⁻¹. We calculate the predicted characteristic finger velocity by combining the buoyancy and mass fluxes to obtain dimensions of velocity, as well as assuming that fingers have a circular cross section (Hoyal et al., 1999b; Carazzo and Jellinek, 2012) V f = g ′ 2 5 ( π V p δ 2 4 ) 1 5 , (13) where V p is given by Eq. 4 and g ′ is calculated using Eq. 12 and a maximum concentration of upper fluid in fingers of 7% (Figure 8). The characteristic length of the instability δ is measured experimentally and estimated from volume conservation to be half the finger width, which is the length scale controlling finger dynamics (Hoyal et al., 1999b). The calculated evolution of finger velocity with initial particle concentration is in good agreement with the observed finger speed (R ² = 0.84).

The evolution of the maximum number of fingers per unit width as a function of both particle concentration and time has been investigated (Figure 9). The number of fingers increases with particle concentration, with 0.25 fingers per cm formed on average at C u = 1 g l⁻¹ and 0.6 fingers per cm at C u = 10 g l⁻¹. This new set of experiments agrees very well with the previous finding for the number of fingers (Scollo et al., 2017), while also expanding the previously-investigated concentration range (Figure 9A). We observe the occurrence of fingers from a few seconds after barrier removal. The maximum number of fingers is achieved 10 to 20 s after barrier removal and this timing weakly depends on the particle concentration (Figure 9B). For high particle concentrations, the number of fingers increases at a faster rate, reaching its maximum value quicker than for low particle concentrations. For example, whilst the number of fingers increases from 0 to a maximum of 0.53 fingers per cm in 10 s for C u = 6 g l⁻¹, it takes 20 s to achieve the maximum number of fingers of 0.27 fingers per cm for C u = 2 g l⁻¹. The number of fingers per unit length gradually decreases after reaching its maximum value (Scollo et al., 2017) and attains a nearly constant value of approximately 0.2–0.3 cm⁻¹, regardless of the initial particle concentration. This value is possibly related to the progressive depletion of particles from the upper layer and to the merging of fingers.

As expected for such settling-driven gravitational instabilities (Hoyal et al., 1999b; Carazzo and Jellinek, 2012), all characteristic length scales, which we measure 10 to 20 s after removing the barrier, are related to each other by a proportionality coefficient. We find finger width and spacing to scale as approximately two and three times the thickness of the PBL, respectively ( W ≈ 2 δ and λ ≈ 3 δ ; Figures 10A–C) and that all these length scales decrease with the initial particle concentration (Figures 10D,E). From the definition of the Grashof number (Eq. 1), the PBL thickness is expected to scale as (Hoyal et al., 1999b; Carazzo and Jellinek, 2012) δ = ( G r ν 2 g ′ ) 1 3 . (14) Finger width and spacing should, therefore, scale similarly. Taking G r = G r c = 10 3 provides an estimate of the critical PBL thickness (Turner, 1973; Hoyal et al., 1999b) above which convection can start (i.e., the minimum PBL thickness). This assumption, however, is found to underestimate the thickness of the PBL measured in experiments. It is therefore probable that the PBL continues to grow beyond this critical thickness. We find that a value G r = G r exp = 10 4 provides a better fit between observations and Eq. 14. The agreement between the measurements and the model is good for the finger width and spacing, whilst Eq. 14 predicts the thickness of the PBL with less accuracy (Figure 10), which is characterised by more scattered data points due to the effect of barrier removal. Moreover, if we fit the measured characteristic length scales to scale with g ′ m , we find m = −0.43 for the PBL thickness, for which the agreement with Eq. 14 is the lowest, whereas for the finger spacing and width, m = −0.38 and −0.35, which are closer to the predicted value of − 1 / 3 .

The particle mass flux through a horizontal plane at depth l located below the upper layer can be calculated in the experiments as F exp = Δ m A Δ t , (15) with Δ m = m t + Δ t − m t the particle mass difference in the region below l measured in a time interval Δ t = 0.1 s, and A = 7.5 × 30.1 cm² the horizontal cross sectional area of the experimental domain. At any time, the mass of particles below l can be calculated by integrating the particle concentration profile (Figure 6) from the bottom of the tank z = b   ( z = 0 corresponds to the height of the initial density interface) to z = l ; m t = A ∫ l b C ( z ) d z . This estimation of the mass does not account for the mass of particles accumulated at the bottom of the experimental tank and therefore only provides a lower limit on the particle mass flux.

Predicting the evolution of the particle mass flux in the presence of fingers is complex because a complete description requires: i) information on the 3D evolution of the finger velocity and size with time and depth, ii) a good knowledge of the duration necessary to produce fingers, and iii) a description of the particle distribution and velocity field inside fingers. It is nonetheless possible to assess the time at which the particle mass flux is expected to increase (i.e., the time at which first particles reach a depth l ). Assuming that particles are coupled with the flow within fingers that form instantaneously at the start of the experiment, first particles can reach a depth l after a time l / V f , with V f the vertical velocity of fingers given by Eq. 13. In contrast, first particles settling individually from the base of the upper layer at their Stokes velocity (Eq. 4) are expected to reach l after a time l / V p .

Figure 11A shows the temporal evolution of F e x p measured at a depth l = 20 cm for C u = 1, 3 and 5 g l⁻¹. Generally, experiments with high particle concentrations result in larger particle mass fluxes that start to increase quicker than in experiments at lower particle concentrations. For all particle concentrations, F e x p increases progressively from when fingers first arrive, before reaching a more stable plateau. The particle mass flux starts to increase close to the time l / V f (15–35 s), which is the expected arrival time for fingers, and much earlier than the expected arrival time of individual particles (143 s). F e x p continues to increase as more fingers reach z = l and fingers become bigger through merging and entrainment of ambient fluid. The latter behavior where F e x p fluctuates around a plateau that occurs once the number of fingers has stabilised (Figure 9B) and overturning has begun in the lower layer, homogenising the particle concentration. To compare different experiments, we average F e x p for 15 s after the plateau is attained and plot this value as a function of C u (Figure 11B). We find that the particle mass flux increases with the initial particle concentration and that a power-law with a fitted exponent of 1.2 can be used to describe the relation between F e x p and C u . In particular, we note that this exponent is independent of l and varies within the interval 1.2 ± 0.1. Theoretically, we can expect F e x p ∝ V f C u . From Eqs 13, 14, the finger velocity scales with C u 4 / 15 so we therefore expect the particle mass flux to be proportional to C u 19 / 15 . The value of the exponent 19/15 ≈ 1.27 is close to the exponent of the power-law used to fit the evolution of F e x p with C u in experiments and there is a good agreement between experimental measurements and the proportionality relation F e x p ∝ 0.14 C u 19 / 15 .

The tendency for settling-driven gravitational instabilities to develop below particle suspensions and to generate fingers has previously been shown to depend on the particle size. Carazzo and Jellinek (2012) compared the instability growth rate with the terminal fall velocity of individual particles and showed that fine ash promotes finger formation, whereas Scollo et al. (2017) experimentally demonstrated that coarse particle (>125 µm) settled individually in water. When particles are sufficiently small, settling-driven gravitational instabilities can form below a particle-laden layer if the time scale of particle delivery to the PBL, given by H c / V p , with H c the thickness of the particle-laden layer, is greater than the time necessary to grow an unstable PBL, given by δ / V p . This condition is independent of the particle size (i.e., V p ). Assuming that the particle concentration is homogeneous in the particle-laden layer, and that its thickness remains constant and uniform, we can therefore assess the potential for instabilities to form by evaluating the dimensionless quantity L ∗ = H c δ . (16)

Combining Eqs. 12, 14, 16, L ∗ can be expressed as L ∗ = H c ( g X p ( ρ p − ρ a ) ρ a G r ν 2 ) 1 3 . (17)

L ∗ >> 1 guarantees that the PBL can reach its characteristic thickness and that fingers can possibly develop. Conversely, L ∗ << 1 means that the particle-laden layer is substantially thinner than the critical PBL thickness, thus preventing finger formation. In our experiments, we consider H c = H 1 , G r =   G r exp = 10 4 and L ∗ consequently ranges from 11.3 to 24.3 in type A experiments, which all produced fingers (Table 2). In order to better constrain the conditions leading to the formation of fingers, we explore a wider range of initial particle concentration and upper layer thickness in type B experiments for which L ∗ varies from 0.8 to 25.

Figure 12 shows a regime diagram for experimental fingers as a function of upper layer thickness and initial particle concentration. Whilst fingers are observed in all the experiments performed at L ∗ ≥ 5, experiments with L ∗ ≤ 2 are characterised by individual particle settling and the absence of fingers (Supplementary Video S5). The domain defined by 2 < L ∗ < 5 corresponds to a transition regime where fingers are either present or absent, depending mainly on the particle concentration, with higher particle concentrations favoring finger formation and their development, since they are more numerous and descending faster. We note that, as expected for settling-driven gravitational instabilities, high particle concentrations and thick upper layers favor the formation of fingers (Jacobs et al., 2015).

The combination of PLIF and particle imaging techniques in separate reproducible experiments provides an approximation of the bulk mixture density evolution with time. We report the formation of an unstable region in all type A experiments, which is more pronounced for high particle concentrations (Figure 6). From experimental observations of the evolution of the particle concentration and of the fluid and bulk mixture density, we summarise the mechanism leading to the onset of settling-driven gravitational instabilities as:

1. The initial density configuration is gravitationally stable, with a light particle-laden layer emplaced above a denser layer;

This mechanism is consistent with theoretical scenarios describing the formation and the destabilisation of a dense particle-laden layer due to particle settling as illustrated in Figure 1 (i.e., the growing of the “nose region”; Burns and Meiburg, 2012; Burns and Meiburg, 2015; Yu et al., 2013; Davarpanah Jazi and Wells, 2020). It is worth mentioning that such a formation mechanism does not require the particle concentration inside the PBL to be greater than the initial particle concentration in the upper layer. The increase in bulk density inside the PBL is due to particles settling across the interface into the upper part of the dense lower layer that therefore becomes heavier than the fluid below. This is in agreement with our experiments where we do not observe an increase of the particle concentration inside the PBL within the uncertainty of our measurements (Figure 6; Table 3). The fact that we do not observe an increase of the particle concentration across the interface is due to the lower layer being at most 1% denser than the upper layer. Thus, in contrast with configurations involving large density differences (Carey, 1997; Manville and Wilson, 2004), particles are not significantly slowed at the interface by this density contrast and do not accumulate. Similarly, the density gradient created at the base of the upper layer (Figure 6A) because of barrier removal only causes negligible variations in the particle settling velocity and does not result in an increase of the particle concentration as reported by Blanchette and Bush (2005) for configurations with significant settling velocity variations.

Furthermore, our experiments suggest that fingers do not inject significant amounts of upper layer fluid into the lower layer and have an average concentration of upper fluid < 3% (Figure 7). This observation indicates that settling-driven gravitational instabilities only weakly affect the transport of fluid phases from the upper to the lower layer, although we observe discrete fluid leaks with concentration of upper fluid up to 7% injected through fingers propagating from the base of the upper layer (Table 3). The small upper fluid concentrations inside fingers suggest that their density, and hence their velocity, depends only weakly on the upper layer fluid density and is mainly controlled by the density of the lower layer fluid and the particle concentration inside fingers.

As in previous studies, we find that, along with particle size, particle concentration exerts a primary control on particle sedimentation and the formation of fingers (Del Bello et al., 2017; Scollo et al., 2017). Our results are in very good agreement with previous measurements of the number of fingers (Figure 9; Scollo et al., 2017) and with theoretical estimates of the finger velocity (Figure 8) and characteristic lengths (Figure 10; Hoyal et al., 1999b; Carazzo and Jellinek, 2012). This increases the confidence in these theoretical formulations of finger geometry and speed on an extended particle concentration range. However, we underline that the scaling relationships of Eqs. 13, 14 are not complete solutions describing the complex flow mechanisms taking place in settling-driven gravitational instabilities. They are based on simplifying assumptions regarding the choice of a single length scale characterising finger dynamics. Moreover, they do not account for the complex particle-fluid and particle-particle interactions happening inside fingers and can therefore give only first order estimates of the size and velocity of both experimental and natural fingers.

We find that the critical thickness determined with G r c underestimates the observed PBL thickness (Figure 10D) and we report that the dimensions of the instabilities formed at different particle concentrations are consistent with a value of G r e x p = 10⁴. This can be interpreted as evidence that the PBL continues to grow even after it reaches its critical thickness ( G r = G r c ). Alternatively, this may also indicate that the effective PBL thickness is affected by barrier removal that initially modifies the flux of particles crossing the density interface by creating small scale eddies that quickly entrain particles below the interface. Even if we assess the dimensions and velocity of the fingers when the effects of the initial perturbation dissipate, this experimental issue can cause differences between the measured PBL thickness and the critical thickness which would develop in an ideal case. Hence, to predict the effective instability dimensions, one should first evaluate G r from observations of the instabilities’ characteristic lengths before using Eq. 14 to extrapolate the dimensions of the instability at different particle concentrations (i.e., different values of the reduced gravity).

We find the convective scaling law of Hoyal et al. (1999b), that works well for experiments, to greatly underestimate the dimensions of volcanic ash fingers by several orders of magnitude. For example, Manzella et al. (2015) estimated the width of ash fingers at Eyjafjallajökull 2010 to be about 170 m from visual observations. Based on observations of the fingers downward velocity, they also estimated the fine ash volume fraction to be in the range of 1 × 10⁻⁶ to 4 × 10⁻⁶ and the PBL bulk density to be approximately 1.31 kg m⁻³. Using the values of the air kinematic viscosity υ = 3 × 10⁻⁵ m² s⁻¹, air density ρ a = 1.30 kg m⁻³ and the Grashof number associated with experiments G r e x p = 10⁴, Eq. 14 only yields a finger width of a few centimetres when applied to the volcanic cloud of Eyjafjallajökull. This suggests that the characterisation of the PBL in nature is highly uncertain. Therefore, this law may not be appropriate for volcanic clouds. In addition to the challenges associated with the measurements of parameters such as the fine ash concentration in volcanic clouds and the PBL, a variety of mechanisms that are not described by our experiments or by those of Hoyal et al. (1999b) and that can form fingers can also explain these discrepancies. In fact, the PBL can purely grow by particle settling, as showed in Figures 1A,B for our experimental configuration, and is most likely to do so in the absence of other vertical fluid motions. Evidence for this mechanisms is perhaps suggested by lidar measurements during the 1991 Mount Redoubt eruption (Hobbs et al., 1991) that showed ash fingers forming at the flat base of the volcanic cloud about 150 km from the vent (“ash veils” in their Plate 1). However, other processes present in volcanic ash clouds can also affect the sedimentation of fine ash. These include jets from overshoot regions, internal waves, overturning motions inside the cloud, and wind-driven stirring (e.g., Carazzo and Jellinek, 2012; Freret-Lorgeril et al., 2020). In addition, Mammatus clouds, i.e., formations characterised by the development of lobes on the cloud base (Schultz et al., 2006), have also been observed at the base of volcanic clouds, e.g., Mt. St. Helens, USA, 1980 (Durant et al., 2009; Carazzo and Jellinek, 2012). However, even if they share visual and dynamic similarities with settling-driven gravitational instabilities, it remains unclear if there is any link between the two phenomena.

Another source of discrepancy with the experiments is related to the fact that, unlike volcanic clouds, experiments are performed in a confined and quiescent environment, without shear at the base of the particle suspension. Although our experiments can address the case where the plume and fingers are advected at wind speed (Scollo et al., 2017), the effect of shear on the base of volcanic clouds needs to be considered in future work in order to describe more adequately the formation and the size of ash fingers. The presence of shear at the base of the cloud can affect the formation of ash fingers for two main reasons: i) shear can inhibit the formation of a PBL by producing eddies at the base the cloud that impede the sedimentation of particles across the interface and ii) it can create Kelvin-Helmoltz instabilities that will interact with settling-driven gravitational instabilities and possibly control the spatial distribution, dimensions and timing of ash fingers. Moreover, our experiments do not explore the effect of internal cloud dynamics on settling-driven gravitational instabilities. Whilst our experiments are in a transitional regime, volcanic ash clouds and ash fingers are fully turbulent and can contain vertical fluctuations that possibly affect the supply rate of particles to the PBL. We expect that more complete observations of ash fingers, combined with dedicated numerical simulations, will provide a better characterisation of ash finger sizes.

Based on a new dimensionless number L ∗ (Eqs. 16, 17), we assess the possibility for the PBL to grow to its characteristic thickness and therefore constrain the conditions favoring the development of fingers. We can distinguish between three regimes in experiments: i) fingers invariably form for L ∗ ≥ 5 whereas, ii) they never form for L ∗ ≤ 2 and iii) are either present or absent for 2 < L ∗ < 5 which corresponds to a transition regime. In general, our analysis indicates that fingers are more likely to form below thick particle-laden layers associated with large particle concentrations, which correspond to L ∗ ≥ 5. Calculating L ∗ in volcanic ash clouds is more complex and requires a good characterisation of the PBL thickness and of the eruptive parameters associated with volcanic eruptions. For instance, measurements of the fine ash concentration in volcanic clouds remain uncertain, as in situ characterisation of volcanic ash is often associated with the very edge of volcanic clouds due to flight restrictions of aircraft (Weber et al., 2012; Eliasson et al., 2014; Fu et al., 2015; Eliasson et al., 2016). Because of these caveats, we are not able to reliably quantify the value of L ∗ for natural volcanic ash clouds.

However, both the cloud thickness and the particle concentrations are related to the MER (Wilson et al., 1978; Sparks, 1986), suggesting that eruptions with high MERs are more prone to develop ash fingers (Scollo et al., 2017). Whilst L ∗ depends on both the particle concentration and the thickness of the upper layer (i.e., volcanic cloud), other parameters also affect the occurrence of settling-driven gravitational instabilities and the production of fingers. For example, particle size has been shown to exert a major control on the tendency to form fingers, with small particles more likely to settle collectively than coarse particles (Carazzo and Jellinek, 2012; Scollo et al., 2017). Hence, although L ∗ provides a good characterisation of the particle concentrations and cloud thicknesses that are associated with fingers, it is not sufficient to obtain a complete characterisation of the conditions favoring the generation of settling-driven gravitational instabilities. L ∗ needs to be considered along with other parameters, such as the particle size and the cloud velocity, in order to assess the potential of volcanic ash clouds to produce fingers.

As with all experimental studies, a number of limitations must be considered when interpreting the data.

First, because of the additional light scattering effect of particles, we only measure the fluid density field in PLIF experiments for low particle concentrations (1 and 2 g l⁻¹), with relatively large uncertainties (±0.8 kg m⁻³ on average below the initial density interface). Future combined PLIF and particle imaging will allow simultaneous measurement of the particle spatial distribution and fluid density (Borg et al., 2001; Dossmann et al., 2016), thus accounting for the effect of particles on fluid density measurements. The combination of these two techniques would therefore contribute to reduce the uncertainty associated with fluid density measurements and provide estimates of the fluid phase properties for particle concentrations >2 g l⁻¹.

Second, measurements of the PBL thickness are initially affected by a disturbance to the density interface generated when removing the barrier. We therefore measure the PBL thickness after it reaches a constant value once the initial perturbation disappears (Supplementary Figure S5). The initial perturbation could be reduced by using a thinner composite barrier (Dalziel, 1993; Lawrie and Dalziel, 2011) and a motorised barrier removal system (Davies Wykes and Dalziel, 2014) in order to obtain a better description of the initial growth of the PBL.

Third, the scaling of the finger size with the Grashof number works well when applied to experiments but not when applied to volcanic ash clouds. As described in Nature and Size of Volcanic ash Fingers, differences between the experimental configuration and volcanic clouds can explain this discrepancy, in addition to the challenges associated with the measurement of relevant parameters in natural eruption. Experiments involving particle-laden currents instead of an immobile particle suspension will contribute to better understand this discrepancy and the scaling of natural ash fingers, along with efforts in numerical simulations of the processes. In fact, the scaling analysis in Scaling of Experiments notably revealed that volcanic clouds are associated with a wider range of dimensionless numbers than our experiments and that R e and G r can be orders of magnitude greater in nature (Table 1). This shows that volcanic clouds are much more turbulent than the analogue particle suspensions in the experiments, that are emplaced in a quiescent environment. Contrarily to the experiments, we can therefore expect turbulent motions inside the clouds to have a more pronounced effect on ash sedimentation, in particular by modulating the flux of particles entering the PBL, as described in Nature and Size of Volcanic ash Fingers. The discrepancy between Grashof numbers associated with the experimental and natural configuration is related to the difficulty of measuring the PBL thickness below volcanic clouds. We expect that a more complete characterisation of the parameters associated with volcanic clouds will contribute to better assess the size of the PBL and evaluate the relevance of the scaling presented in our experimental study.

Finally, laboratory experiments, unlike volcanic plumes, exist in a confined space affected by boundary conditions and the formation of return flows (Roche and Carazzo, 2019). In our experiments, boundary conditions have an impact on convection in the lower layer, which has only a limited effect on measurements of finger characteristic velocity and dimensions, especially in the early part of the experiments, and a negligible influence on finger formation. Despite these uncertainties however, we can infer useful and rigorous interpretations of the results.