Rheology of three-phase suspensions determined via dam-break experiments

Three-phase suspensions, of liquid that suspends dispersed solid particles and gas bubbles, are common in both natural and industrial settings. Their rheology is poorly constrained, particularly for high total suspended fractions ($\gtrsim 0.5$). We use a dam-break consistometer to characterize the rheology of suspensions of (Newtonian) corn syrup, plastic particles, and CO$_2$ bubbles. The study is motivated by a desire to understand the rheology of magma and lava. Our experiments are scaled to the volcanic system: they are conducted in the non-Brownian, non-inertial regime; bubble capillary number is varied across unity; and bubble and particle fractions are $0\leq\phi_{gas}\leq0.82$ and $0\leq\phi_{solid}\leq0.37$ respectively. We measure flow-front velocity and invert for a Herschel--Bulkley rheology model as a function of $\phi_{gas}$, $\phi_{solid}$, and the capillary number. We find a stronger increase in relative viscosity with increasing $\phi_{gas}$ in the low to intermediate capillary number regime than predicted by existing theory, and find both shear-thinning and shear-thickening effects, depending on the capillary number. We apply our model to the existing community code for lava flow emplacement, PyFLOWGO, and predict increased viscosity and decreased velocity compared with current rheological models, suggesting existing models may not adequately account for the role of bubbles in stiffening lavas.

Three-phase suspensions, of liquid that suspends dispersed solid particles and gas bubbles, are common in both natural and industrial settings. Their rheology is poorly constrained, particularly for high total suspended fractions ( 0.5). We use a dam-break consistometer to characterize the rheology of suspensions of (Newtonian) corn syrup, plastic particles, and CO 2 bubbles. The study is motivated by a desire to understand the rheology of magma and lava. Our experiments are scaled to the volcanic system: they are conducted in the non-Brownian, non-inertial regime; bubble capillary number is varied across unity; and bubble and particle fractions are 0 ≤ φgas ≤ 0.82 and 0 ≤ φ solid ≤ 0.37 respectively. We measure flow-front velocity and invert for a Herschel-Bulkley rheology model as a function of φgas, φ solid , and the capillary number. We find a stronger increase in relative viscosity with increasing φgas in the low to intermediate capillary number regime than predicted by existing theory, and find both shear-thinning and shear-thickening effects, depending on the capillary number. We apply our model to the existing community code for lava flow emplacement, PyFLOWGO, and predict increased viscosity and decreased velocity compared with current rheological models, suggesting existing models may not adequately account for the role of bubbles in stiffening lavas. dies show that the apparent viscosity of particle suspensions increases with particle volume fraction, and intermediate and concentrated suspensions show non-Newtonian behavior with finite yield stresses and strain-rate dependence [2,3,4]. Bubble suspensions exhibit an increase or decrease in apparent viscosity relative to the liquid phase. Concentrated bubble suspensions also show non-Newtonian behavior, including non-zero yield stresses and shear-thinning or shearthickening behavior [5,6,7,8,9,10,11,12,13,14,15]. A key parameter determining the rheology of bubble suspensions is the capillary number: where µ is the liquid viscosity, a is the undeformed bubble radius,γ is the shear strain rate, and Γ is the surface tension. The capillary number reflects the ratio of viscous stresses that deform bubbles and capillary stresses that restore them [7,13,14] When Ca 1, viscous stress dominates, bubbles deform easily and decrease relative viscosity; when Ca 1, the surface tension resisting bubble deformation dominates, bubbles remain nearly spherical and increase relative viscosity. The rheology of three-phase suspensions (liquid suspending particles and bubbles) has received much less detailed investigation, but is known to be a complex function of the suspended phase fractions, micro-textural properties such the size and shape distributions of particles and bubbles, and the conditions of shear [10,16,17]. The only systematic experimental study, which was limited to φ solid 0.5, φgas 0.3 and low capillary number [17], found that the rheology was well described by a simple convolution of existing two-phase rheology models in which the liquid and bubbles were treated as an effective medium that suspended the particles. We expect that interactions between bubbles and particles will become increasingly important at higher bubble and particle fractions such that a simple convolution of two-phase models is insufficient.
Magma (and lava, which is its subaerial counterpart) is a natural three-phase suspension, composed of a molten silicate liquid (melt) that suspends a variable fraction of solid particles (crystals) and gas bubbles. The suspending melt is Newtonian over a wide range of strain rates, and has a viscosity that can vary over orders of magnitude [18,19]. Bubble and crystal volume fractions can both range from 0 to 1, and typically change, along with melt viscosity, during transport through the Earth's crust and over its surface. The rheology of magma exerts a firstorder control on the dynamics of its eruption, and the emplacement of subsequent lava flows [20,21]. Whilst the rheology of natural magmas has been measured directly [22,23,24,25], such measurements are challenging to perform and interpret, and are not well suited to systematic investigation of parameter space. As a result, it has been common for rheology of magmatic suspensions to be investigated via analogue experiments [4,6,13,17,26]. Even with the use of analogues, it has proven challenging to perform rheometry on samples with φgas 0.5. This is because samples with high bubble fraction are difficult to prepare, particularly when the sample also contains particles, and because they are prone to breakdown when loaded into a conventional rheometer. We conduct experiments in which bubbles are grown in situ via a chemical reaction, and analysed using a dam-break consistometer. This approach circumvents both problems, allowing us to investigate samples with bubble fractions (0 ≤ φgas ≤ 0.8) and particle fractions (0 ≤ φ solid ≤ 0.37) that span the most relevant ranges for natural magma and lava.

Theoretical background (a) Scaling
We scale our analogue suspensions to basaltic magma, which typically has a pure-melt viscosity in the range 10 2 to 10 4 Pas (10 3 − 10 5 Poise) at eruption temperature [19], and surface tension between 0.05 and 0.3 N/m [27,28,29]. Bubbles that nucleate deep within the volcanic plumbing system are likely to be small (radius ∼ 10 −5 to 10 −3 m) and experience very low shear rates (∼ 10 −10 to 10 −5 1/s), resulting in capillary numbers Ca 1 [30,31]. As magmas decompress on their way to the surface, bubbles grow and coalesce to radii ∼ 10 −5 − 10 −2 m or larger [32,33,34]. However, decimeter sized bubbles often segregate from the flow, making a suspension rheology approach no longer appropriate. Magma in basaltic dikes experiences strain rates up to ∼ 10 1 to 10 2 1/s [35] and capillary number therefore spans from Ca 1 to Ca 1. At the surface, lava flows, particularly fast-moving flows (velocities ∼1 m/s) may experience strain rates of 10 −3 to 10 2 1/s, and both high and low capillary numbers [8]. The viscosity of silicate melts is sufficiently high that viscous stresses at the particle scale dominate over inertial and Brownian stresses [4]. We therefore conduct experiments that span from low to high capillary number, whilst remaining in the non-inertial, non-Brownian regimes.
(b) Theoretical framework for rheology of three-phase suspensions To our knowledge, the only theoretically-grounded model for the rheology of three-phase suspensions is that of Phan-Thien and Pham [16]. They present a model for the effective Newtonian viscosity, η, of a three-phase suspension of the form: where ηr,s is the relative viscosity of the three-phase suspension (i.e. suspension viscosity normalized by the liquid viscosity), and ηr,p and η r,b are the relative viscosities of the twophase components (particle-liquid, and bubble-liquid suspensions respectively). Importantly, this relation treats one two-phase system as an effective viscous medium in which the third phase is suspended. Truby et al. [17] validated this theoretical form against experimental data, using existing models for the viscosity of the two-phases suspensions: following, respectively, Mueller et al. [4] and Llewellin and Manga [7], where Einstein exponents B solid = 2 and Bgas = 1, and φm is a maximum packing fraction of solid particles. Both Phan-Thien and Pham [16] and Truby et al. [17] note that the definition of φgas and φ solid depends on which two-phase suspension is chosen as the effective medium; in this study we choose the particle suspension such that φ solid = V solid /(V solid + V liquid ) and φgas = Vgas/(Vgas + V solid + V liquid ), where V denotes the volume of the subscript phase. Truby et al. [17] extend the model of Phan-Thien and Pham [16] to allow for non-Newtonian effects by equating the suspension viscosity, η, with the consistency, K, in the Herschel-Bulkley model [36]: where τ is the shear stress, τy is the yield stress, and n is the flow index (n > 1 implies shearthickening, n < 1 shear-thinning). We adopt this approach in the analysis below.

Methods (a) Experimental setup
Our experiments use a Bostwick (dam-break) consistometer, in which fluid is released from a reservoir and into a confined rectangular channel (Fig. 1   To the syrup, we add rigid plastic particles sieved to a mesh size of 60-80 (Sourced from Precise Finishing Inc Stock No. GP-PP60/80, Polyplus / Type II). Particles have near-neutral buoyancy (ρsyrup = 1395 kg/m 3 , ρ plastic = 1370 kg/m 3 ), are angular and slightly elongate (aspect ratio 1.5 to 2), and have a dominantly unimodal size distribution (long axis ≈ 200 to 400 µm). Particles are stirred into the suspension by hand which leads to minimal air entrapment as shown in the microphotograph in Supplementary Fig. S2. Volume fraction of solids in the suspension is determined prior to the introduction of gas by mass of syrup and particles (uncertainty of 1-2 volume %). When multiple experiments with the same particle content were performed, a large volume of syrup and particles was prepared, but gas was introduced to each sample individually.
Gas is introduced via a chemical reaction between baking soda (NaHCO 3 ) and citric acid (C 6 H 8 O 7 ) which produces water vapor and CO 2 gas. Reactants are added with a ∼3:4 ratio by mass, yielding an empirical relation of roughly 0.5 g baking soda and 0.7 g citric acid per 100 g syrup, per 10% desired gas volume fraction. The reaction begins immediately and the sample is allowed to rise for ∼20 min prior to its release into the confined channel. Samples with φgas < 0.6 are risen in a separate container, then transferred to the consistometer reservoir; samples with φgas > 0.6 are more delicate, and are risen in the reservoir to avoid foam disruption during transfer. Final gas volume fraction is calculated from the final mixture density. Mass of the suspension is measured when transferring to the channel and the volume is calculated using the internal channel dimensions and the height of the mixture within the reservoir as measured through the channel walls. This method leads to ∼3 to 5% relative uncertainty on the volume and resultant ∼3 to 5% relative uncertainty on the density and an absolute uncertainty of ∼0.02 to 0.04 on φgas.
The water vapor that forms during the reaction is insufficient to lower the syrup viscosity appreciably. In the experiments with the highest concentration of reactants (9 g total per 100 g of syrup), that water could, in principle, result in a dilution from 77% to 76.6% sugar if all the water dissolved, which would change the viscosity from 7.6 Pas (76 Poise) to 6.5 Pas (65 Poise) at 23 • C. However, a test suspension, from which the evolved gas was allowed to escape overnight showed no additional dissolved water when analysed using the Brix refractometer.
The flow begins with the removal of the dam, which takes less than one second. The dam is typically coated with ∼2 to 3 mm of syrup, when removed, effectively shortening the length of the fluid in the reservoir. However, in this geometry, propagation is insensitive to the length of the reservoir and this is below the spatial resolution of our forward model. Flow progression along the channel is recorded by time-lapse photographs from the side at a rate of 1 frame-per-second (FPS), and from above using a video camera (30 FPS). Experiments last between 15 and 180 s, with most 30 to 90 s in duration. We process the top-view videos to extract the flow front position over time. In the tracking algorithm, the video, rotated to align the flow direction left to right, is processed by manually identifying the starting position, scale, and center line in each video, and setting an appropriate threshold on the brightness of either the red, blue or green channel of the image to distinguish syrup from the background. In some cases, the time-distance profile of the flow front was smoothed to remove artifacts caused by the flash from the side-view camera. The flow front position is compared to a forward model of dam-break flow to invert for the rheological parameters.

(b) Extracting rheological parameters from experiments
We follow the derivations of Liu and Mei [39] and Balmforth et al. [40] for time evolution of flow thickness, h, of a Herschel-Bulkley fluid flowing down a slope: where x is the along-channel distance, ρ is the fluid density, g is the acceleration due to gravity, and θ is the channel slope. This forward model is valid for non-inertial laminar flow of a gravitationally settling fluid experiencing no-slip conditions at the base of the flow and no-stress at the free surface. In the inertial regime, propagation of the flow front would be expected to scale with the shallow water speed √ gH. where L is the reservoir length and H the reservoir depth, we expect inertia would be dominant for t 10 −5 − 10 −2 s, which is much shorter than the duration of our experiments. We confirm that all experiments are in the laminar regime by verifying that the Reynolds number, Re, is below the cutoff for the onset of turbulence at Re ∼ 1000. We calculate Re for a free surface flow of a Herschel-Bulkley fluid in a rectangular channel as Re = [41] whereū is the mean velocity, D is the hydrodynamic radius, D = 4HW 2H+W , for channel width W = 15 cm, and channel depth H ∼ h. We find Reynolds numbers 10, well below the cutoff. We confirm no-slip conditions via imaging through the base of the channel, which shows that bubbles are well coupled to the base ( Supplementary Fig. S3). We neglect the effect of surface tension between the syrup and air, which becomes important for h < 2 Γ ρg ≈ 0.5 cm [42]. Accordingly, we restrict our experiments to have a starting reservoir volume H > 4cm. Finally, we require that the suspension be able to deform under its own weight, which requires that the Bingham number, B < 1, where B ≡ τyL ρgH 2 [40]. We prepared a suspension of φ solid = 0.46, and find that for a flow height of 2.5 cm, it does not deform under it's own weight, suggesting that yield stress must exceed 42 Pa.
We solve equation 3.1a numerically using a finite-difference scheme on a centered 3-point stencil which is second-order accurate in space and we use two-step Runge-Kutte for second order accuracy in time. Initial conditions are fluid height h = H where x ≤ L (within the reservoir), and h = 0 elsewhere. Boundary conditions are no-slip on the bottom, no-stress at the free surface, ∂h ∂x |x=x L = 0 (no inflow on low slope), and h|x=x R = 0 where x R is chosen to be sufficiently far from the final position of the flow front such that flow does not reach the far boundary (domain length is 20% longer than the final experimental flow length). The domain is discretized in space evenly with 52 grid points and with a constant time step of ∆t = 6.67 × 10 −4 s. The forward model assumes that temperature, rheology, and density fluid of the fluid are constant in time and space. We allow slope to be non-zero in the model solution to account for sensitivity of flow advance rate to slope angles below our uncertainty for direct measurement (≤ 0.5 • ). We extract the position of the flow front over time from the model by identifying where h ≈ 0, below a threshold of 1 × 10 −4 to account for numerical diffusion, smoothed linearly between grid points in x.
We use the evolution of the flow front over time in the experiments as a constraining observation in our inversion for rheological parameters. We invert for rheological parameters (K, τy, n) and the slope (θ) using an Ensemble Kalman Filter (EnKF) approach [43,44]. This is a probabilistic approach in which the forward model is instantiated many times with varying parameters that are modified iteratively to minimize the misfit between the model and data (typically ∼2000 data points). All samples are propagated together between iterations, which allows faster convergence than a random walk. We use Gaussian priors centered about a visuallydetermined close fit for the rheological parameters, and about zero for the slope. Standard deviations were 10%, 10%, 5%, and 1 • for the initial distributions of K, τy, n, and θ, respectively. EnKF simulations are initialized with 300 samples and run for 5 iterations; convergence is typically achieved within 4 iterations, insensitive to small changes in the initial parameters ( Supplementary Fig. S5). In other Earth-science applications, EnKF is often used to update probabilities by assimilating more information through time, for example in volcano monitoring [43,44]. We do not utilize this capability of EnKF; all information about the flow progression is included from the beginning. We use the Python implementation of the EnKF algorithm provided on GitHub by geoyanzhan3 1 [43,44]. We report the numerically-determined maximum likelihood value, and uncertainty given by the 5% and 95% quantiles. We find uncertainty values of 20% in K, 20% in τy, and 10% in n.  Table S2.

Results
We present results from 35 dam-break experiments spanning 0 ≤ φ solid ≤ 0.37 and 0 ≤ φgas ≤ 0.82 ( Fig. 2 and Table S2). We find that, with increasing volume fraction of solids, relative consistency (Kr = K/µ) increases, yield stress increases, and flow index decreases (more shearthinning), consistent with previous studies [1]. Increasing gas fraction in all cases increases relative consistency and yield stress. Gas-bearing experiments show both shear-thinning (n < 1) and shear-thickening (n > 1) behavior, with increasing n correlated with increasing Ca. We use our measurements to parameterize the model framework presented in Section 2b, and find empirical expressions for Kr, τy, and n as functions of phase volume fractions and (where relevant) the maximum capillary number. We calculate maximum capillary number using equation 1.1 with the maximum observed bubble size from each experiment (0.5 to 8.2 mm) and estimate the maximum strain rate from experiment videos. Best fitting parameter values for these empirical models are determined using the L-BFGS-B algorithm for limited-memory, quasi-Newton, bound-constrained optimization [45] implemented in Python; uncertainty is estimated from the approximated Hessian matrix.
We find the following empirical relationships: where φm=0.56±0.20, B solid =2.74±1.56 and Bgas = 1.98±0.09, C 1 = 80.0±10.9, C 2 = 1.98±0.23, the critical volume fraction for the onset of appreciable yield stress φc,τ y = 0.35±0.01, C 3 = 0.70±0.25, C 4 =0.55±0.31, and the critical volume fraction for the onset of shear rate dependence φc,n = 0.39±0.12. Results from each experiment are plotted in Fig. 2 and compared to contours for the best-fitting model plotted on the same color scale. Misfit between the experiments and model are given in Table S2.

Discussion
(a) Comparison with previous studies (i) Relative Consistency, K r Our expression for consistency (Eq. 4.1a) has the functional form presented in Eq. 2.1, such that it can be decomposed into separate functions for the two-phase components, facilitating comparison with previous studies (Fig. 3A and 3D, equations and parameters in Table S4). For particle suspensions (φgas = 0), we find a Kr(φ solid ) dependence that has the same functional form as that of Krieger and Dougherty [3] and Roscoe [46]. We do not have experiments for φ solid > 0.37 where relative consistency rapidly increases and a transition to another functional form, such as that in Costa et al. [47], may be appropriate. For bubble suspensions (φ solid = 0) we find a Kr(φgas) dependence that is stronger than most previous studies, with bubbles having a stronger stiffening effect on the suspensions. Furthermore, we find that in our experiments, bubbles always increase the relative consistency, regardless of capillary number, in contrast with previous studies [6,13]. We note that our consistency model (Eq. 4.1a) is optimized against the full three-phase dataset, hence three-phase interactions are likely to play a role in these discrepancies with models that are fitted against two-phase data. For instance, breakup of bubbles between interacting particles, and the formation of liquid-film bridges between particles may both amplify the effect of the bubble fraction on suspension rheology, and cause the rheology to be controlled by the smallest bubbles, which are in the low Ca regime, despite observations of deformation in the largest bubbles (Fig.  S2), which imply Ca 1 locally.
(ii) Yield Stress, τ y We select an empirical exponential-type relationship for yield stress (Eq. 4.1b), which overlaps with previous two-phase studies ( Fig. 3B and 3E). We observe a much stronger effect of particles compared to gas, but onset at around the same volume fraction. Yield stresses we observe in threephase suspensions are higher than would be predicted by a linear superposition of the two-phase components, indicating again that interactions among the phases are important. At intermediate volume fractions we find moderately good agreement with Hoover et al. [49] and Mueller et al. [4]. The gas-liquid suspension experiments compare favorably to results by Princen and Kiss [11] and Rouyer et al. [12], although there is limited overlap between regions of calibration.
(iii) Flow Index, n We fit flow index data using a piecewise continuous linear relationship after Castruccio et al. [26] and Truby et al. [17]. For particle suspensions, and for bubble suspensions at low Ca, this formulation yields a shear-thinning behavior, which is consistent with previous relationships (Fig.  3C and 3F). With increasing Ca and gas fraction, it transitions to shear-thickening, a behavior reported in previous studies [10].

(b) Data limitations
We use the dam-break consistometer setup to measure the rheology of three-phase suspensions in order to avoid bubble breakdown which is common in conventional rheometers. This setup allows us to attain high bubble volume fractions, but requires numerical inversion and additional considerations for direct comparison to rheometry data collected via other methods. The behavior of Herschel-Bulkley fluids in confined channels and on unconfined slopes has been investigated theoretically and experimentally [50], primarily using kaolin clay suspensions and Carbopol which show variable agreement with measurements collected in conventional rheometers. Most studies that use channel flow to investigate rheology use the approach to a final runout state [51], the height profile [52,53,54], the internal velocity [55,56,57,58], or the runout distance with time    [40,52,59,60]. The internal velocity and height profiles observed in experiments using kaolin clay typically compare favorably to theory, except during the initial slumping stage, in which inertia is dominant, and near the flow front where significant curvature requires more sophisticated modeling [53,58]. In contrast, experiments using Carbopol show systematic offsets compared to conventional rheometers, possibly a result of basal slip, three-dimensional effects such as levee formation [55,56], or a scale effect introduced in the narrow gap of conventional rheometers compared to a relatively large microstructural length scale in Carbopol [58]. Our results use the runout distance with time to measure rheology which shows good agreement with conventional rheometry for experiments on corn syrup, xanthan gum, and kaolin suspensions [40]. We expect our suspensions to behave similarly to two-phase kaolin particle suspensions, and further verify that basal slip and levee formation do not occur in our experiments. Given findings in previous work on two-phase suspensions in dam-break experiments and the good agreement between our experiments and existing particle-bearing suspensions, we do not expect systematic offsets in the rheometric parameters determined in our study compared to what would be measured in a conventional rheometer. Our forward model is 1D (depth-integrated) and cannot capture any effects perpendicular to the channel direction. Visual inspection of the experiments suggests that there is not pronounced cross-channel flow. Our channel width, at 15 cm, is larger that that used by Balmforth et al. [40] who find a 10 cm channel outperforms a 5 cm channel. We estimate the boundary layer width  Table S2. The misfits between these experiments and the model do not show systematic differences nor are the misfits different in magnitude compared to experiments with low predicted boundary layer thicknesses.
The uncertainties on the values for Kr, τy, and n we derive from inversion of the experimental data are 25%, 31%, and 14%, respectively, and the mismatch with models for those parameters (Eqs 4.1a -4.1c) are 35%, 106%, and 10%. Our experiments show good repeatability at φgas ≤ 0.6. If we compare experiments for which φ solid and φgas are within ±5%, and de-trended Kr, τy, and n values using their respective models, we find ranges of 38%, 149%, and 20% in those parameters, which are of the same order as uncertainties and mismatch. For φgas > 0.6, the data show greater variability with ranges of 68%, 154%, and 73%, respectively. Data variability may be related to experimental sources of uncertainty, such as small errors in flow density or height of fluid in the reservoir to which the forward model is very sensitive, or microstructural differences in the particle and bubble populations. Yield stress measurements at all bubble fractions show scatter between experiments and misfit to the model higher than the uncertainty predicted by EnKF model fitting, which may come from real variation among samples rather than trade-offs between fitting parameters.
We see consistent deviations between the experimental data and the rheology model in the range of φgas > 0.6, where the model consistently over-predicts both the relative consistency and flow index, perhaps due to trade-offs in the inversion of the experimental data. To test the importance of trade-offs in fitting parameters, we find the Pearson correlation coefficients between the EnKF instances after the final propagation (best fit for the fitting parameters) for each experiment. A value of 1 indicates a perfect positive correlation, -1 a perfect negative correlation, and 0 no correlation.The EnKF analyses reveal a moderate correlation between the slope angle and consistency, with an average correlation coefficient between EnKF samples for each experiment of C=0.63. We also find weak correlations between slope angle and yield stress (C=0.25) and slope angle and flow index (C=-0.44); a weak positive correlation between consistency and flow index in experiments with φgas > 0.6 (C=0.27); a moderate negative correlation when φgas < 0.2 (C=-0.40); and no apparent correlation at intermediate gas volume fractions 0.2 < φgas < 0.6 (C=-0.01). We do not observe correlations in the parameter inversions between yield stress and consistency (C=-0.06) or yield stress and flow index (C=0.11). For an example of the EnKF inversion results and match to experimental data see Supplementary Fig. S5.
Our results consider strain rates up to 10 1 1/s and maximum Ca between ≈ 10 −2 and 10 1 . We find that bubbles lead to a higher effective viscosity than predicted by the Cross model [1], but our experiments do not cover a large-enough range of Ca to comment on asymptotic behavior that may occur at much higher or lower Ca. Given the strong strain-rate dependence of the Herschel-Bulkley formulation, this model may not be appropriate for fluids experiencing strain rates far outside the calibration region ( 10 1/s) as may occur in natural settings, including during conduit ascent in explosive eruptions or in fast-moving lava flows.
Additional uncertainty exists for the estimate of Ca, which is difficult to define for a polydisperse bubble population whose exact bubble sizes and deformation cannot be directly observed. Estimates for the bubble size distribution are available in the supplementary material. Bubble deformation is observed for some of the largest bubbles experiencing the strongest shear, and thus experiencing the highest Ca conditions, imaged through the side walls of the channel (Supplementary Fig. S2). Such observations are limited to the exterior of the flow, which is necessarily dominated by edge effects. Better characterization of bubble sizes and deformations in the flow interior is an avenue for future work.

(c) Application to lava dynamics
We demonstrate the implications of our new three-phase suspension rheology model with realistic parameters for lava flows through implementation to the channelized lava flow model PyFLOWGO [62,63]. PyFLOWGO is a 1D lava flow advance model for the velocity and final runout distance of a cooling and crystallizing lava. The physical model uses a Bingham-plastic rheology of the form τ = τy + η ef fγ . We add our experimentally constrained rheology model to several built-in parameterizations for the yield stress and effective viscosity. We calculate τy using Eq. 3b and the effective viscosity term η ef f = Krµγ n−1 using Eqs. 3a and 3c. We compare our model with the commonly used two-phase modified Einstein-Roscoe relationship [46], and three-phase model of Phan-Thien and Pham [16].
We take the flow simulation parameters from the well-documented LSF1 lava flow from Mount Etna, Italy in 2001. We initiate the simulated lava with a low solid volume fraction of solid = 0.1. Solid fraction then increases, with further crystallization during emplacement, to φ solid = 0.15 when using our rheological relations, and φ solid = 0.3 for the modified Einstein-Roscoe and Phan-Thien and Pham [16] relationships. Gas fraction begins and remains at φgas = 0.65. The simulated flow follows the path of steepest descent in the pre-eruptive topography. The initial flux is 20 m 3 /s. We use a temperature-dependent melt viscosity using the VFT equation (a = -4.827, b = 5997, c = -330.3) [64, 65,66]. We tune the model parameters of initial channel width (set to 40 m) and temperature (set to 1223.15 K) to obtain the observed flow height (10 m) and mean velocity (10 −3 to 4 × 10 −2 m/s) when using our new rheology model. The width and initial temperature are held constant across all simulations. To predict the effective viscosity, we use a reference strain rate of 10 −3 s −1 for this lava flow, estimated from observed velocities and flow thickness. We approximate the shear rate dependence using the capillary number calculated for a 1 cm radius bubble and a surface tension of air in basalt of 0.37 N/m [29]. As the flow cools, the melt viscosity, crystal content, and effective bulk viscosity are updated every step.
Our model predicts a higher effective viscosity compared to that of Phan-Thien and Pham [16] or Roscoe [46] by more than an order of magnitude, and yields a mean velocity that is ∼3 times slower (Fig. 4). Additionally, our yield strength model diverges from that based on Dragoni [67] ∼80 m down flow, when the crystal fraction increases enough to support significant bubble-crystal interactions. Although shear-thinning behavior is not included in PyFLOWGO, we compare effective viscosities predicted by our model assuming a Bingham-plastic rheology (n = 1) and with a dynamic value for n in response to changing crystal fraction (while maintaining a constant shear rate). We highlight that for the large bubble volume fraction of some natural lavas, n deviates from one (initially n ≈ 0.82) and can result in a factor of two difference in effective viscosity, and a 15% difference in runout distance . This effect may be more pronounced in modeling that includes a physical model for shear-rate dependence in the flow.

Conclusion
We use a dam-break setup to measure the rheology of three-phase suspensions, scaled to apply to magma and lava. Through the use of a chemical reaction between baking soda and citric acid, we are able to create corn syrup suspensions with bubble volume fractions up to 0.82 and particle volume fractions up to 0.37. Strain rates in our experiments range up to 10 1 1/s and maximum Ca between ≈ 10 −2 and 10 1 , relevant to basaltic lava flows. We invert for a Herschel-Bulkley rheology using a depth-integrated 1D forward model [39] and the probabilistic Ensemble Kalman Filter approach.
Based on our results, we develop a suite of models for the Herschel-Bulkeley parameters, calibrated across the whole dataset, which are presented in equations 4.1a-4.1c. Combined with equation 2.3, this constitutes a rheological model for three phase suspensions. Our results show good agreement with existing literature for particle-bearing suspensions with a strongly nonlinear increase in viscosity with increasing particle fraction, development of yield stresses at particle fractions φ solid 0.35, and development of shear-thinning behavior at particle fractions  Pham [16] (three-phase). We also highlight the contrast between assuming n=1 (Bingham-plastic rheology) and the Herschel-Bulkley rheology (n≈0.8) that results in a decrease in effective viscosity by a factor of 2 at low crystal fractions.
φ solid 0.39. Our findings indicate that suspended gas leads to a non-linear increase in viscosity with increasing bubble fraction, larger than predicted by previous work, and development of yield stresses and shear rate-dependent behavior at bubble fractions, or combined bubble and particle fractions, similar to those observed for particle-only suspensions. We find that shear-rate dependence of bubble-bearing two-and three-phase suspensions are shear thinning in the low capillary number regime and become shear thickening with increasing capillary number. This study contributes an explicit incorporation of three-phase interactions into a novel model for yield stress, and a parameterization for the inclusion of capillary number in flow index. We highlight the implications of our model for magma/lava dynamics through application to the lava flow emplacement model PyFLOWGO, and show higher viscosities and lower predicted velocities for surface lava flows as a result of the stiffening behavior of concentrated bubble suspensions at low and intermediate capillary numbers. Our findings underscore the need to incorporate three-phase rheology into simulations of magma or lava flow used to address fundamental science questions as well as hazard assessment and mitigation. We characterize the bubble size distributions of our experiments using optical imaging at two different spatial scales. First, we consider the smallest size fraction of bubbles in undeformed syrup suspensions using optical microscopy which allows us to resolve bubbles between 20 µm and 2 mm in diameter. Second, we characterize the size distribution and shear-induced elongation of bubbles 1.5 mm that are visible without magnification through the walls of the clear channel.
The smallest bubble size fractions are determined from photomicrographs analyzed using a MATLAB circle-finding algorithm and FOAMS (Shea et al. 2010). Results show a power-law distribution of bubble number density, shown in Fig. S1 for an undeformed sample with only gas suspended in liquid (no solids). In the presence of particles, heterogeneous nucleation is observed with small bubbles ( 50µm in diameter) nucleating and growing on the uneven surfaces of the particles. This population shows greater spatial heterogeneity, obvious in panel B of Fig. S1 in which a region with fewer particles shows fewer and larger bubbles, while an adjacent area shows a greater density of small bubbles. The high density of bubbles on the particle surfaces introduce significant errors to the size distribution estimate for three-phase suspensions, making quantitative comparison challenging. We might additionally expect the population size to vary in response to shear, but deformation and coalescence at this size fraction are hindered by the relatively high surface tension.
To characterize the larger bubbles visible through the side walls of the channel, we again utilize FOAMS (Shea et al. 2010) to calculate bubble volume, aspect ratio, and elongation direction distributions. We compare the bubbles found in the undeformed reservoir with the flow front after 8 min to show a change in the elongation direction of large bubbles (Fig. S2). The reservoir and the flow front show similar size distributions with a logarithmic decrease in bubble number density with increasing size at a comparable rate to the smaller size fraction (measured on two different sample suspensions). Both populations show slight elongations (mean aspect ratio ≈ 0.6), but the bubbles in the reservoir show a preferential vertical orientation as bubbles rise due to buoyancy, while the bubbles at the flow front are aligned with the flow direction, an approximately 90 • rotation.
We also add a photograph in Fig. S3 taken from the underside of the channel looking up at the flow. It shows bubbles in contact with the bottom of the channel, without the formation of a segregated bubble-free layer. Fig. S4 shows the relative error between experiments and the empirical model, normalized by the experimental results. We also include an example experiment inversion including the posterior distributions of Kr,τy,n,and θ which show strong covariance between the slope and relative consistency, a good fit between observations and the maximum likelihood model, and convergence of the EnKF inversion after 3 iterations. Finally, we report each relationship shown in Fig. 3. Fig S5 shows [2] 1 + Bφ B = 1.5 − 5 Roscoe [46] (1 − bφ) −B b = 1 − 1.35