2 of 19 volume fraction of solidified bubbles, in natural flows varies in wide range, e.g., 0.02-0.88 at Kīlauea and 0.04-0.71 at Etna (Cashman et al., 1994;Chevrel et al., 2018;Gaonac'h et al., 1996;Halverson et al., 2021;Harris & Rowland, 2015;Herd & Pinkerton, 1997;Patrick et al., 2019). Higher vesicularity is measured close to the vent and decreases with flow distance (Cashman et al., 1994). Bubble coalescence in a highly vesicular paho'eho'e lava flow generates various void structures (Walker, 1989). Crystal content in lavas at Kīlauea is 0.05-0.35 (Gansecki et al., 2019;Robert et al., 2014). In addition, magmas and lavas can have strain-rate dependent behavior which can have important consequences on flow dynamics, e.g., by causing the transition from paho'eho'e to 'a'a (Sehlke et al., 2014).
As a fluid dynamics problem, the morphology and velocity of viscous flow on a slope have been investigated extensively theoretically and experimentally (e.g., Craster & Matar, 2009;Didden & Maxworthy, 1982;Huppert, 1982b;Mériaux et al., 2022). An analytical solution for a viscous flow on a sloping plane predicts that the flow length increases with time according to a powerlaw relationship (Lister, 1992), whose power index depends on the fluid's rheological properties (Balmforth et al., 2007).
Bubbles in magma can increase or decrease the effective viscosity depending on whether bubbles can deform under a certain shear rate (Llewellin & Manga, 2005;Manga et al., 1998). Large-scale lava flow experiments using molten basalt have shown that bubbles increase the effective viscosity (Soldati et al., 2020), while field measurements of viscosity at high shear rate suggest the decrease of effective viscosity (Chevrel et al., 2018(Chevrel et al., , 2019)). We discuss this in detail in Sections 2 and 6.
Bubbles and crystals can also sometimes cause a yield strength, a minimum stress required to allow deformation. If the stress caused by the weight of the fluid layer does not exceed the yield strength, it will not flow. Here, the yield stress is mostly due to the microstructure of the suspensions, which can evolve with time. In a fluid layer slightly overcoming the yield strength, a slow downflow changes the microstructure in suspension, decreases the viscosity, and then accelerates the gravity flow (Coussot et al., 2002). Later injection of the low particle fraction fluid causes acceleration of the flow (Castruccio et al., 2014). Yield stress fluids show shear localization even in homogeneous stress fields (Møller et al., 2008). The yield strength and fluid properties depend not only on the particle volume fraction but also on the particle size distribution (Castruccio et al., 2010).
The shear conditions and phase concentrations, and therefore rheology, of lava flows can vary greatly in space and time during an eruption. For example, the voluminous lava flow eruption in 2018 at the lower East Rift Zone (LERZ) of Kīlauea Volcano, Hawai'i, exhibited rapid flow with velocities up to 15 m s -1 within the confined channel (Dietterich et al., 2021;Patrick et al., 2019). Yet overflowed lava from the channel showed slow and thin patches. Such a significant change of lava behavior only a short distance from the source is different from previously observed lava flows (Biass et al., 2019).
At the flow front of some thin overflows, the lava formed large cavities with a length scale of >1 m inside the closed flow patch, called shelly paho'eho'e (Jones, 1943;Swanson, 1973;Wentworth & Macdonald, 1953). Figure 1 shows a photograph of a solidified shelly paho'eho'e formed during the 2018 eruption at the lower East Rift Zone (LERZ) of Kīlauea Volcano, Hawai'i. The standard explanation for forming such cavities is an emptying out of the lava beneath a solidified crust, similar to a small lava tube. However, this flow patch, and others like it, do not have an outlet through which to discharge the lava (Figure 1). This feature requires an alternative mechanism for formation in which gas accumulates at the toe of the flow. Here, we propose a formation mechanism through a process of flow localization and large bubble concentration at the flow front. We investigate this mechanism, and the behavior of highly vesicular flows in general, through a series of experiments with bubbly lava flow analogs.
Bubbles and particles in fluid suspensions change the effective viscosity of suspensions, relative to that of the liquid phase, η l . We here define the bubble volume fraction ϕ b and the particle volume fraction ϕ p as
where V l , V b , and V p are the liquid, bubble, and particle volumes, respectively. We define the particle volume fraction with respect to the dense rock equivalent (bubble-free) suspension to be consistent with the literature (e.g., Chevrel et al., 2018;Gurioli et al., 2014).
The effective viscosity of a liquid-particle flow increases with the particle volume fraction. One such model expresses the viscosity ratio η r as
where η eff is the effective viscosity of the particle-bearing fluid, ϕ p is the volume fraction of particles, and ϕ max is the maximum possible volume fraction of particles or undeformed bubbles in random close packing, and the exponent is 2.0 ≤ Cϕ max ≤ 2.5 (Mader et al., 2013). The random close packing of spheres of uniform size is ϕ max = 0.64. This value decreases with increasing particle aspect ratio and increases with increasing size variation (Klein et al., 2017;Mueller et al., 2011). An alternative view of the particle size contribution is that the small particles together with the liquid phase form a viscous fluid phase, and the large particles determine the solid fraction in Equation 3 (Castruccio et al., 2010).
The effective viscosity of bubbly fluids depends on the capillary number (Llewellin et al., 2002;Manga & Loewenberg, 2001;Pal, 2003), which is a measure of whether the bubbles will get deformed by the flow. It is defined as the ratio of viscous stresses acting to deform bubbles, and surface tension resisting deformation
where ̇ is strain rate, R is the bubble radius, and Γ is the surface tension. In this equation, the viscosity η can be either the effective viscosity η eff or liquid viscosity η l . When a large bubble is surrounded by a suspension with smaller bubbles and particles, the dominant viscous force for the deformation of the bubbles depends on the effective viscosity η eff . In contrast, the viscous force to deform an isolated bubble should be defined with the liquid viscosity η l . For Ca < 1, bubbles do not deform and behave similarly to solid particles, thus increasing the effective viscosity as is expected by Equation 3 by replacing ϕ p to the volume fraction of bubble ϕ b . In this case, ϕ max in Equation 3 becomes the random close packing of undeformed bubbles, and for monodisperse bubbles ϕ max ∼ ϕ bmax ∼ 0.64. Again, ϕ bmax is sensitive to the bubble size distribution.
For Ca > 1, bubble deformation becomes significant and decreases the effective viscosity. Some theoretical and empirical equations have been suggested for the impact of bubbles on the effective viscosity for different Ca regimes. A numerical estimate of the relative viscosity for an infinite shear rate (high Ca) is close to the Hashin-Shtrikman upper bound (Manga & Loewenberg, 2001) We here emphasize quantitative estimates of effective viscosity based on above models have limitations. Bubbles within a suspension can have significant size variation, increasing the maximum bubble fraction and modulation of effective viscosity for undeformed bubbles (Equation 3). Bubble size variation also alters the criteria of deformation through Ca (Equation 4). Temporally evolving bubble shape and size by coalescence can affect the effective viscosity of bubbly fluids (Takeda et al., 2020). A qualitatively meaningful interpretation of the above scalings is that bubbles can increase the effective viscosity by several orders of magnitude while decreasing it by only 1 order of magnitude (Figure 2).
Bubbles are deformable to be polyhedral so that the bubble fraction can exceed the random close packing of monodisperse spheres (ϕ b > ϕ bmax = 0.64). Such a bubbly fluid has elasticity by surface tension. An empirical relation describes the shear modulus as a function of the bubble radius and fraction (Hohler & Cohen-Addad, 2005)
Here, ∼ indicates approximation in a factor.
For a rhyolitic magma, the effect of surface tension is negligible relative to the shear modulus of the liquid phase, the shear modulus decreases with bubble fraction
where G 0 is the shear modulus of bubble-free suspending fluid and 1 < ζ < 4 is constant (Namiki et al., 2020).
The yield stresses τ y reported for foams by oscillatory measurements are approximately described by the following relation (Hohler & Cohen-Addad, 2005):
(10)
When both bubbles and particles exist in a suspending fluid and bubbles do not deform, the volume fraction of the deformable suspending fluid decreases. Thus, adding bubbles to a dilute particle suspension at the low capillarity number regime increases its viscosity. In contrast, adding bubbles to a concentrated particle suspension, where space of the deformation is limited, decreases its viscosity (Truby et al., 2015). An effect of high-volume fraction of particles on the behavior of bubbles has been experimentally confirmed (Pistone et al., 2016). The three-phase fluid also has elasticity. The shear modulus can be written by the product of bubble effect (Equation 8) and particle effect (Equation 3; Namiki & Tanaka, 2017)
Recall that Equation 3 accounts for reducing the viscously deformable space.
Equation 11 has similar dependence on ϕ p , suggesting a reduction of the elastically deformable space.
Birnbaum et al. ( 2021) developed a set of empirical expressions to describe the rheological parameters of three-phase suspensions with both particles and bubbles that is based on a Herschel-Bulkley formulation (Herschel & Bulkley, 1926)
where τ is the shear stress, τ y is the yield stress, K is the consistency, and n is the flow index (n > 1 implies shear-thickening, n < 1 shear-thinning). Using analog dam-break flow experiments in confined channels, the following expressions were derived:
where ϕ m = 0.56 ± 0.20, B p = 2.74 ± 1.56 and B b = 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.
The common characteristic of the three-phase fluid is that the solid particles can increase viscosity and elasticity significantly, while the effects of bubble is limited (Birnbaum et al., 2021;Namiki & Tanaka, 2017;Truby et al., 2015).
To understand the effect of a high-volume fraction of bubbles and crystals on a lava flow, we conducted a series of experiments of a particle-bearing bubbly lava flow analog, made of corn syrup as the liquid phase, plastic particles as crystal analogs, and CO 2 bubbles. We simulated the lava flow by pouring the lava analog on an inclined acrylic plank (Figure 3). The experimental conditions are summarized in Table 1.
To create our lava analog, we mixed corn syrup and plastic particles in a beaker until the mixture was uniform and then added reactants into the mixture. The chemical reaction caused a volume increase of the mixed suspension. The corn syrup had a temperature and water-content dependent viscosity (Figure S1 in Supporting Information S1) that we varied between 0.35 and 150 Pa s and surface tension of 0.08 N m -1 (Rust & Manga, 2002). We conducted all experiments at the same room temperature so that the temperature-dependence of the syrup viscosity did not impact the results. We present the water-content dependence at the same temperature in Figure S1 in Supporting Information S1). The plastic particles had a density in a range of 1,400-1,500 kg m -3 , were irregularly shaped with an aspect ratio of 1.5 and 0.1 mm in size in averages, and served as analogs to crystals.
We introduced bubbles by a chemical reaction of baking soda and citric acid, which generated CO 2 . This is the same technique and materials used in previous works (Birnbaum et al., 2021;Namiki & Tanaka, 2017). NAMIKI ET AL. 10.1029/2022JB024139 6 of 19
The amount of bubbles could be roughly controlled by adjusting the amount of baking soda and citric acid. However, the gas produced could not be finely controlled because it escapes into the atmosphere during the reaction. The bubble volume was calculated by subtracting the volume of liquid and particles from the total volume of the suspension. We measured the total volume of the three-phase suspension using a graduated beaker. By introducing bubbles, the liquid color became whitish, indicating that the fluid color can be a measure of bubble fraction. The liquid and particle volumes were calculated by their mass and density. The vesicularity range used in our experiment was 0.65-0.78, which simulated the highly vesicular lava flow (0.5-0.82) that erupted during the 2018 LERZ eruption from Ahu'ailā'au (formerly known as Fissure 8) at Kīlauea Volcano, Hawai'i (Patrick et al., 2019; Figure 1). Crystallinity in lavas from Ahu'ailā'au is low (0.05), while that for the early phase is 0.35 (Gansecki et al., 2019). The particle fraction in our experiments simulated the latter case.
After the volume increase stopped and before any significant liquid separation from the foam by gravity occurs, we poured the fluid from the beaker by hand onto a plastic plank, inclined at angles of 25-28° from horizontal. This slope angle was chosen to be consistent with conditions in nature and allow measurements of active flow with the given materials. We tried slopes from 8° to 45° and found that a flow on a steeper slope lengthens in the direction of the downslope, while a flow on a shallow slope acquires a round shape. We also found a technical problem that the flow down a shallow slope of <20° is too slow, and its surface dried during the experiment due to water evaporation. A very steep slope >30° is rare in nature (Biass et al., 2019;Favalli et al., 2010).
The duration of pouring was on the order of 10 s. We tried various methods of introducing the particle-bearing bubbly mixture onto the plank, such as injection through a tube or removal of a dam, and found that controlling the flow rate of the non-Newtonian sticky fluid was quite difficult. In addition, flow in a tube modified the bubble fraction in the fluid before it reached the plank. We thus concluded that this simple method of pouring directly from a beaker is the best way as a first step to simulate short lived unconfined bubbly lava flows.
We observed the analog lava flow on the transparent acrylic plank using video cameras, and quantified the velocity field at the surface on the videos using optical flow velocimetry (Lev et al., 2012;Sun et al., 2010). For the measurement of the velocity field of pure liquid flow, we used a more viscous fluid than that of the reference experiment. This was because to measure the velocity field for a pure liquid, we needed to add tracers on the flow surface. Adding tracers to the rapid flow with low viscosity is challenging. Note that the appearance of two pure liquid flows was similar in terms of no localized flow. In some experiments, we used a laser displacement sensor (Omron ZG2, vertical resolution of 6 μm) to measure the surface topography of the flow. We also used a more viscous liquid for topography measurement. The surface of the bubbly flow with lower viscosity liquid had more large bubbles, which became scattering sources of the laser beam. We calculated the time evolution of the topography at two representative locations from the line sensor output. Characteristics of bubbles can vary with time and space, which cannot be monitored quantitatively. In this series of experiments, we thus try to understand the bubbly flow in a qualitative manner.
Slope Plank length Liquid viscosity Liquid density Total volume ϕ b ϕ p (°) (m) (Pa s) (kg m -3 ) (L) (-) (-) Figure 4a 28 1.2 0.35 1,315 0.44 0.00 0.00 Figures 4b and 6b 28 1.2 0.35 1,315 0.43 0.69 0.00 Figures 4c and 6c 28 1.2 0.35 1,315 0.42 0.69 0.35 Figure 6a 25 0.66 150 1,470 0.50 0.00 0.00 Figure 6d 26 0.66 1.9 1,370 0.38 0.65 0.00 Figure 7 25 0.66 1.0 1,350 0.10 0.78 0.00 Note. Total volume is the poured volume on the plank for the model lava flow. Table 1 Experimental Conditions Journal of Geophysical Research: Solid Earth NAMIKI ET AL. 10.1029/2022JB024139 7 of 19
We can estimate the yield stresses for our bubbly lava analog. According to Equation 10, the yield stress for bubbly fluid depends on the maximum bubble volume fraction. Assuming ϕ bmax = 0.64, our bubbly fluid with ϕ b ≤ 0.78 and bubble size of >10 -4 m has yield stress of <10 Pa. Using Equation 8, the shear modulus is <100 Pa, resulting in that the critical strain for yield becomes <0.01. Similarly, the yield stress for our three-phase suspension is estimated using Equations 13b, 24 Pa, and the shear modulus is 200 Pa (Equation 11), resulting in the critical strain for yielding of 0.1. Note that the characteristics of the three-phase fluid, such as liquid surface tension and particle size, can affect also the yield stress. The fluid used by Birnbaum et al. (2021) to derive Equation 13bwas approximately the same as that in this paper.
Figure 4 summarizes three typical experiments. We used three kinds of fluids; pure liquid, bubbly fluid, and particle-bearing bubbly fluid. The most significant result is that the pure liquid flows rapidly, while the bubbly fluid and the particle-bearing bubbly fluid flow slowly. At the downslope nose of the pure liquid, droplets originating from capillary instability precede the flow (Figure S2 in Supporting Information S1; Huppert, 1982a;Lister, 1992).
To visualize this trend quantitatively, we plot the time evolution of the length of each flow (Figure 5a). Theoretically, the length of flows would be expected to increase with time t in a power-law manner t α -for a Newtonian fluid, while the fluid is being poured at a constant flux, = 7 9 , and after the pouring ends, = 1 3 , as denoted in Figure 5a (Lister, 1992).
The pure liquid, expected to behave as a Newtonian fluid, flows down rapidly and the flow length time evolution follows the power index of 7 9 , as expected for a Newtonian fluid. This result suggests that the flux rate we achieve by pouring is nearly constant, despite pouring by hand. Given that the flow is too rapid, the flow front reaches the bottom end of the plank before the pouring ends, and we could not obtain the power index of length evolution after the pouring ends (Figures 4a and 5a).
In contrast, the flow of the bubbly fluid advanced much more slowly than the pure liquid (Figure 4b). Bubbly fluids are often non-Newtonian, and advance according to different power-law indexes (Balmforth et al., 2007;Birnbaum et al., 2021). In our experiment for bubbly fluids, the power-law index before the end of pouring deviates from the = 7 9 expected for Newtonian fluids. This may suggest the effect of non-Newtonian rheology, but NAMIKI ET AL. 10.1029/2022JB024139 8 of 19 because of short time pouring different from Balmforth et al. (2007), the power index until the end of pouring does not reach a constant. Once the pouring stopped, though, the flow advance rate follows a power-law index of 1 3
, which is consistent with a Newtonian behavior (Figure 5a). The behavior of the particle-bearing bubbly fluid is more complex. The time evolution of the length approximately stops around 100 s, but the fluid resumes to flow down around 500 s at a rate following a power-law index of 1 3 , consistent with Newtonian behavior. In Figure 4c at 17 min, the width of the flow front is narrower than that of the original, indicating that the former entire flow upslope adheres to the plank and only the central part of the flow is mobile. The width of the flow front decreases with time. The low power index before reaching
3 may be due to the rheology effects such as yield stress and shear rate dependence (Balmforth et al., 2007;Birnbaum et al., 2021). The time varying power index may suggest that these rheological parameters are also time dependent.
The maximum widths of these three examples do not differ significantly (Figure 5b). The bubbly fluid shows a slight decrease in measured maximum width with time after the pouring end. This is because the fluid upslope at the maximum width flows down to make the measured width narrower. In addition, some bubbles at the edge of the flow disappear with time. The width of particle-bearing bubbly suspension continuously increases for a long time. This is because the initially thick flow spreads by gravity in both the downslope and cross-slope directions.
The width of flows is also expected to increase over time in a power-law manner, t α . For a Newtonian fluid, while the fluid is being poured at a constant flux, = 1 3 , and after the pouring ends, the growth stops and α = 0 (Lister, 1992). In Figure 5b, the initial phase is faster than = 1 3 , which may originate from the width of the beaker, which is wider than a point source. The later stage of pouring is close to = 1 8 . The power index of 1 8 is derived, assuming that fluid flux is approximately zero at the very beginning of pouring. Maybe our hand moved downslope. This should halt the width evolution at a certain location but not much affect the evolution in length.
In bubbly fluid experiments (Movie S1, Figure 4b), a dark linear structure appears along the center axis at the end of the pouring, which overtakes the other whitish parts downslope. Here, the dark fluid color is due to fewer bubbles, implying that this region is denser than the surrounding whitish region. We here quantify this feature by measuring the flows' surface velocity field (Figures 6a-6c) and flow thickness (Figure 6d). Figures 6b and 6c is Figures 6a-6c show that the flow front has higher velocity than upslope. The pure liquid flow, as a Newtonian fluid, has a limited velocity variation in the cross-slope direction (Figure 6a), while the bubbly fluid shows strongly localized flow (Figure 6b). This region corresponds to the dark linear structure, indicating that the dense region with fewer bubbles flows faster than the surrounding. The particle-bearing bubbly flow shows that the flow is concentrated in the narrow flow front, which corresponds to the region that resumed motion after coming to a near halt (Figure 6c).
Motivated by the importance of thickness on flow velocity, we then measured the flow thickness using a laser displacement sensor (Figure 6d). The liquid viscosity in Figure 6d is higher than that used in Figure 4b, and bubbles in Figures 4b and 6d maintain spherical shapes, suggesting that Ca < 1. The inset shows a photo of the measured flow, in which darker spots appear at the middle of the flow. These darker spots, together with the whitish region between them, rapidly flow down. Initially, the flow front is upslope of the sensor. After 40 s, the flow front reaches the location of the sensor (Figure 6d). As expected, the flow front is relatively thick, and the height decreases with time as the flow is advected downslope. So, the x-axis time in Figure 6d can be considered as the distance from the flow front. The center of the flow (red curve) is thicker than the side (blue curve) at the flow front, but later the side becomes thicker than the center.
When we made the bubbly fluid in the beaker, the liquid was segregated by gravity and formed a lower liquid layer with dark color (Figure 3, red arrow). This process may also be important in the flow. We then observed the bubbly flow from the bottom of the plank (Figure 7). The bottom view shows a darker color, suggesting that more liquid accumulates at the bottom of the flow. The bottom view also shows that, as the experiment progresses, the darker region concentrates at the flow front. The side view indicates the flow front is thicker than the upslope region. The top view shows an increased concentration of bubbles at the flow front. Journal of Geophysical Research: Solid Earth NAMIKI ET AL. 10.1029/2022JB024139 10 of 19
Our experiments show that the bubbly flow is slower than the pure liquid flow for the same liquid viscosity. Recall from Section 2 and Figure 2, bubbles can increase or decrease the effective viscosity. Because, in our experiments, bubbles only slowed the flow velocity, we infer that the bubbles did not substantially deform.
To confirm this hypothesis, we estimate the capillary number for our bubbly fluid flow. The vertically averaged mean flow velocity in the layer u for a Newtonian fluid on a slope angle of θ is calculated by (Turcotte & Schubert, 2014)
where ρ = ρ l (1 -ϕ b ) is the effective density of the fluid, g = 9.8 m s -2 is the gravitational acceleration, h is the thickness of the flow. Equation 14gives a simple order of magnitude estimate of velocity, which in reality varies throughout the flow depth and with any spatial or temporal variations in the flow. Indeed, our experiments show that the flow thickness and bubble distribution vary locally and temporally (Figures 6 and 7). In a natural lava flow, the slope angle also differs locally. We regard Equation 14 as applicable for the local average in a relatively homogeneous region.
We take the average shear rate ̇ in the definition of Ca (Equation 4) to be u/h. Substituting Equation 14 and this definition of ̇ into Equation 4, the viscosity term disappears and we obtain
In this equation, the liquid density, surface tension, and slope angle do not vary significantly. That is, the flow thickness and bubble size are the dominant factors determining the capillary number. In our experiments, the flow thickness varies temporally and locally and is on the order of 10 -2 m (Figure 6d). The bubble sizes also vary in the range of 10 -4 -10 -2 m. The minimum size of 10 -4 m is measured under microscope (Figure 3), and the maximum is measured from videos of analog lava flow, taking the time variation into account (Figures 4 and 7).
We plot Equation 15in Figure 8 and find that the possible bubble radius and flow thickness range in our experiments are located in the field of Ca < 1, indicating that the bubbles do not deform and effective viscosity increases. This estimate is consistent with our observation that the bubbles always slow the flow velocity. Journal of Geophysical Research: Solid Earth NAMIKI ET AL. 10.1029/2022JB024139 11 of 19
We here note that Equation 15 assumes that a large bubble is surrounded by bubbly fluid. When ϕ b is small or bubble size is uniform, liquid viscosity may be more appropriate to be a substitute in Equation 4. Given that most bubbles did not deform in our experiments Equation 15becomes
The last term decreases with increasing ϕ b to reduce Ca. Thus, estimated Ca in Figure 8 is the maximum estimate. Again, our experiments are located in the field of Ca < 1 irrespective of the surrounding fluid properties.
In Section 4.2, we quantified the characteristics of flow localization. In the bubbly flow, the line-like dark region has a thinner thickness and flows faster than the surrounding region (Figures 6b and 6d). To understand the source of this structure, we revisit the initial condition. When we made the bubbly fluid in the beaker, the lower part of the beaker became relatively liquidrich by gravitational segregation (Figure 3, red arrow). As a result, when we poured the bubbly suspension on the plank from the beaker, the fluid that came out later had fewer bubbles. Since the driving force of the downslope flow is gravity, this later fluid parcel, which has fewer bubbles and thus a higher bulk density and lower effective viscosity, flows down faster than the surrounding fluid and becomes thinner. These characteristics are consistent with those observed in Figures 6b and 6d. We thus interpret this region to have a lower bubble fraction, originating from liquid separation in the beaker. Indeed, Movie S1 shows that the line-like dark structure appears at the end of the pouring. Our interpretation is consistent with the previous work in which later injected less viscous fluid accelerates the flow (Castruccio et al., 2014).
Such an artificial bubble distribution does not occur in nature. However, the temporal fluctuation of bubble fraction in an erupted lava has been observed during the LERZ eruption (Patrick et al., 2019). The above result indicates the importance of the heterogeneous bubble distribution in the gravitationally driven downflow. Next, we consider the effect of the vertical localization of bubbles, which may occur in natural systems.
In experiments with particle-bearing suspensions (Figure 4c), the flow resumed after almost stopping. It is known that suspension fluids with yield stress begin to flow abruptly once a threshold stress is overcome (Coussot et al., 2002). This is because the rheology depends on the internal microstructure and deformation fluidizes the jammed suspension to accelerate the flow.
The estimated yield stress for our particle-bearing bubbly fluid is 24 Pa (Equation 13b), which is close to the stress caused by the fluidal static pressure acting on the fluid layer ρ l (1 -ϕ b )gh sin θ ∼ 20 Pa with a thickness of h ∼ 10 -2 m. This agreement explains the observation that the particle-bearing bubbly flow almost stops (Figure 5a, red curve, 100 s). The fluid was initially thick and generated the stress exceeding the yield stress ρ l (1 -ϕ b )gh sin θ > τ y . With time, the fluid becomes thinner by flowing, and the stress acting on the flow falls below the yield stress, and then the flow stops. From our visual observation, it is not plausible that the critical strain of 0.1 for yielding occurs within the three-phase suspension to resume the downflow after the stopping. However, the flow resumed.
In contrast, the yield stress of the bubbly fluid of <10 Pa is much smaller than the fluidal static pressure, which is consistent with the observation that the bubbly flow did not stop (Figure 4b). Irrespective of differences in yield stress, both flow lengths had the same power index of = 1 3 close to that typical of a Newtonian fluid at the very end of observations (Figure 5a). 15). The red and blue lines indicate Ca = 1 for a lava flow and model lava flow in our experiments, respectively. The orange and blue area indicates a possible range of real lava flow and model lava flow, respectively. We here assume ρ l = 2,700 kg m -3 , θ = 25° (solid line), θ = 1° (dotted line), Γ = 0.4 N m -1 for a lava flow, and ρ l = 1,400 kg m -3 , θ = 25°, Γ = 0.08 N m -1 for experiments. For lava flow, the range of thickness and slope is estimated by the eruptions at Kīlauea (Biass et al., 2019;Dietterich et al., 2021;Patrick et al., 2019) and Mt. Etna (Favalli et al., 2010). The bubble size is based on the field observations (Chevrel et al., 2018;deGraffenried et al., 2021;Shea et al., 2010; Figure 1). For experiments, the range of the flow thickness is estimated based on the side view (Figure 7) and the measurement by the laser displacement sensor (Figure 6). The bubble radius range is estimated by the analog lava flow video and microscope measurements (Figures 3, 4, and 7).
We interpret this phenomenon as the result of vertical bubble localization developed with time. Less dense bubbles are buoyant in a bubbly suspension, and bubbles and liquid separate by gravity with time. Indeed, Figure 7 shows that the liquid in the bubbly flow concentrates at the bottom of the flow. Bubbles floated by buoyancy, and the liquid segregated to accumulate at the bottom.
Summarizing the above discussion, we interpret our experiments as shown in Figure 9. Liquid concentrates at the bottom of the flow front above which the bubbles float. The observation that our experiments approach the Newtonian scaling with time suggests that the bubble-poor regions at the base and center of the flow control the velocity of flow front propagation. Based on this idea, we estimate the time scale of vertical bubble separation.
The density contrast between liquid and bubble causes bubble floating and liquid drainage. We here focus on the two-phase system of bubble and liquid for simplicity. When the bubble fraction is lower than the random close packing ϕ b < ϕ max , the evolution of liquid-rich layer (dashed line in Figure 9) is produced by the floating of bubbles.
An isolated bubble in a liquid layer ascends under the balance of buoyancy and viscous drag which is known as Stokes velocity (Turcotte & Schubert, 2014)
In a suspension, the ascending of individual bubbles will be hindered by the presence of neighboring bubbles. This effect is empirically tested with settling particles (Davis & Birdsell, 1988). Applying this relation for bubbles, the hindered rise velocity becomes
Actual growth rate of the liquid-rich layer may be intermediate of Equations 17 and 18. This is because, in the experiment on gravity-driven settling of packed particles into a liquid layer, the growth rate of the liquid-rich layer is intermediate between velocities estimated by these two equations (Shibano et al., 2012).
Another way to estimate the liquid separation is using Darcy velocity, in which liquid phase flows down in the suspension
where we assumed that pressure difference is originated from the density difference of mixture and liquid phase. At a high fraction of bubbles, the bubble fraction varies with depth, as denoted in Figure 9, and the drainage in such foam is complicated (Cantat et al., 2010). Coarsening and coalescence of bubbles occur, which vary the bubble size with time and depth (Maurdev et al., 2006;Sato & Sumita, 2007). Compaction induced by the bubble deformation is also expected (Bercovici & Ricard, 2003;Huber & Parmigiani, 2018;McKenzie, 1984). These processes are essential to generate the large cavity at the top of the bubbly flow, but the film rupture process of the viscous surfactant-free liquid is not well established (Ohashi et al., 2022), and we here use Equations 17-19 for an order estimate.
Next, we estimate the time scale to develop the liquid-rich bottom layer. In the experiment Figure 4b, ϕ b = 0.69, which is higher than ϕ bmax , but considering the vertical variation of the bubble fraction, it is appropriate to use Equation 18 for a minimum estimate. Using the parameters for the experiment in Figure 4b, assuming a bubble radius of 1 mm and a bubble fraction at the base of the bubble-rich layer ϕ b ∼ 0.64, which is appropriate for a wet foam (Cantat et al., 2010), we obtain v h ∼ 0.07 mm s -1 . Dividing the bubble size by the hindered rise velocity, we estimate that a liquid layer with a thickness of the bubble size develops in ∼10 s. Figure 5a indicates that the power index of the bubbly flow becomes close to 1 3 , a value close to Newtonian fluid, after ∼10 s elapsed. This may be because deformation is concentrated in the bottom Newtonian fluid layer.
The liquid-rich layer can explain the limited effect of bubbles on slowing the flow velocity. The bubble fraction of ϕ b = 0.69 for the experiment in Figure 4b can increase the viscosity 100 times relative to pure liquid (Figure 2). The density has also decreased to about 30% of the pure liquid. Together, these should lead to slowing down the flow <1/100. However, the velocity reduction is ∼1/10 (Figure 4b). The non-Newtonian characteristic might be a source of this discrepancy, but the power index of 1 3 suggests that the liquid-rich layer behaves as a Newtonian fluid (Figure 5a). One possible interpretation is that the bottom liquid-rich layer helps to lubricate the bubbly flow, and the increase of effective viscosity is less than the above estimate. In our particle-bearing bubbly fluid, the particle density is close to that of the liquid phase, and therefore the particles do not segregate as fast as the bubbles. We also note that the volume fraction of particles is less than the close packing (i.e., 0.35 < ϕ max = 0.64), which is too small to generate large yield stress. Thus, the vertical separation of bubbles could cause the flow to resume (Figure 4c). If the particles are denser than the liquid phase or its volume fraction is close to ϕ max = 0.64, the bottom of the bubble-poor layer concentrates with particles, which increases the effective viscosity.
The above discussion provides estimates of the orders of magnitude rather than a quantitative analysis. This is because our bubbly fluid has size distributions that affect the effective viscosity, and the quantitative estimate and control of bubble properties are quite difficult.
Figure 8 demonstrates that we can estimate Ca for individual bubbles in a flow going down a slope only from the fluid thickness and bubble size. This is because the viscosity cancels out, and other parameters are relatively constant. In our experiments, Ca < 1, and bubble deformation is not significant (Figure 8). In this section, we use Figure 8 to explore the range of Ca expected for natural lava flows.
The possible range for real lava flows (orange region) straddles both areas of Ca < 1 and Ca > 1. In particular, low viscosity basaltic flows close to the source region, in confined channels, over steep slopes, or that are several meters thick, can experience rapid velocities and high shear rates and high Ca. Indeed, field viscosity measurements of paho'eho'e magma using rotational viscometer indicate that deforming bubble decreases the effective viscosity from 700 to 380 Pa s, where the sample has vesicularity of 0.5, bubble size of 0.34 mm, and the crystallinity of few vol.% (Chevrel et al., 2018). This measurement was conducted under high shear rate of >1 s -1 , which may be a rare condition within a natural lava flow. We here note that the larger bubbles deform more easily (Equation 4), and larger bubbles are abundant in a real lava flow. In such conditions, even a flow with a small shear rate may cause viscosity decrease as a consequence of bubble deformation.
We here provide a possible example in the field for viscosity decreases and increases associated with changes in Ca. Figures 10a and 10b show an overspill from the main lava flow from Ahu'ailā'au (Fissure 8) vent. The overspill leaves the main lava flow from the edge of the levee at 1.9 km downstream from the vent, and forms fingers after ∼100 m of flow. The velocity of the overspill in this narrow channel is >0.5 m s -1 but in the fingers, it is ≤0.05 m s -1 . Here, the lava flow velocity depends on the flow thickness, slope angle, and viscosity (Equation 14). The daytime view of the solidified overspill (Figure 10c), suggests that the flow thickness and slope do not change with location significantly. Cooling and increasing crystallinity are major controlling factors of lava viscosity in general (Harris et al., 2009;Lipman & Banks, 1987), yet here the overspill emplacement time (∼200 s) is too short to allow cooling to play a significant role in slowing the flow.
One possible interpretation is that the change in Ca, and the associated increase in viscosity, are responsible for the slowing of the overspill flow. The melt viscosity is estimated as 64 Pa s from the viscosity measurements (Soldati et al., 2021). Bubbles and crystals increase the effective viscosity to 4.9 times the liquid viscosity (Gansecki et al., 2019), resulting in the three-phase bulk viscosity becoming 290 Pa s. With bubble radius of 20 mm, Γ = 0.4 N m -1 , shear rate of 0.1 s -1 as a ratio of flow velocity 0.5 m s -1 to the half of channel width 5 m, the estimated Ca ∼ 1. If the depth of the channel is shallower than 5 m, which is plausible, the shear rate associated with the vertical velocity profile is higher than 0.1 s -1 , so smaller bubbles can deform. Large bubbles (R > 0.02 m) exceed Ca > 1 and can deform, which decreases the effective viscosity and accelerates the flow. Such a rapid flow increases Ca, deforms bubbles, and enhances coalescence (Okumura et al., 2006), which may facilitate the formation of giant bubbles. The larger bubbles are more buoyant and concentrate at the top of the flow front. In contrast, smaller bubbles R < 0.02 m are below Ca < 1 and increase the effective viscosity. Here, smaller bubbles may be left in the lava because of slow ascending velocity. Consequently, the effective viscosity of the lava flow increases at the flow front. If the shear rate becomes too high, the flow may disrupt the giant bubbles and generate smaller ones (Stein & Spera, 1992), contributing to reducing Ca suddenly.
Different from our experiments, fresh less viscous lava is supplied from the upstream, which could generate a local viscosity contrast in a lava flow to cause viscous fingering (Bischofberger et al., 2014;Chen et al., 2017;Homsy, 1987;Nittmann et al., 1985;Oppenheimer et al., 2015). The outer shell of each finger may be eventually solidified and the ascending giant bubbles accumulate beneath it. When the bubbles are concentrated, the chance for bubbles to come into contact with each other increases. Assuming that the bubble coalescence is governed by percolation, the increase of local bubble fraction is expected to enhance the coalescence of bubbles to make larger bubbles (Gaonac'h et al., 2003;Giachetti et al., 2019;Rust & Cashman, 2011). Larger bubbles generated by coalescence accelerate gravitational drainage in the bubble film (Nguyen et al., 2013). Drainage makes the bubble film thinner than critical for rupture. Thus, the closely located bubbles eventually coalesce to make a huge cavity at the end of the flow beneath the solidified crust.
As shown in Figure 1, a large cavity with a length scale of >1 m exists without an outlet, suggesting that initially the flow contained voluminous bubbles and then bubbles concentrated at the flow front. The measured vesicularity of the lava flow erupted during the 2018 LERZ eruption Kīlauea Volcano, Hawai'i is 0.5 ≤ ϕ b ≤ 0.82 (Patrick et al., 2019). If bubble separation occurs in a patch of lava with a volume of 1 m 3 and an initial vesicularity of 0.82 becomes finally 0.5, a cavity with a volume of 0.64 m 3 emerges. Thus, the vertical bubble separation likely generated a shelly paho'eho'e.
The lava flow cools from the surface, so bubble separation needs to occur before the solidification. Lava flows erupt at high temperatures, such as the basaltic lava at Hawaii (>1,000 °C). The large temperature difference with the ambient environment cools the lava flow. Heat loss from the lava surface is caused by radiation and convection of surrounding air (Griffiths, 2000).
In Figure 11, we estimate the time scale to develop a solidified crust by the thermal diffusion, where a hot lava is contacting an atmosphere with a constant temperature
where L c is the thickness of the cooled lava crust to a given temperature, t is time, T is the glass transition or a certain temperature, T a is atmospheric temperature, T e is eruptive temperature. We here estimate the thermal diffusivity by the method used in Namiki et al. (2021)
where k m = 1.22 W m -1 K -1 and k a = 0.166 W m -1 K -1 are thermal conductivity of bubble-free magma and steam, respectively, and C p = 2,300 J K -1 kg -1 is heat capacity of melt. We here neglect the effect of bubbles on heat capacity because its effect based on the mass fraction is quite small. The measured vesicularity for the LERZ eruption is 0.5-0.82 (Patrick et al., 2019). For this vesicularity range, the range of thermal diffusivity is narrow 2.2 × 10 -7 < κ < 3.2 × 10 -7 m 2 s -1 , so we used an average value of 2.6 × 10 -7 m 2 s -1 . 10.1029/2022JB024139 16 of 19
The surface temperature of the lava flow is higher than the atmospheric temperature (Fink & Griffiths, 1990). The hot lava heats the air above it and drives air convection, which cools the lava, and then the surface temperature adjusts to the atmosphere (Moitra et al., 2018;Stroberg et al., 2010). Our estimate assumes the constant atmospheric temperature of the lava surface in Equation 20, which implies neglecting the time for surface cooling, and results in a maximum estimate of the solidified crust growth. Convection may occur inside the lava flow, enhancing the heat transport from the inside lava flow to the atmosphere. However, a viscous boundary layer develops at the surface of the lava, in which heat transportation is governed by thermal conduction (Griffiths et al., 2003).
The black solid curve in Figure 11 shows the glass transition temperature of 700 °C (Ryan & Sammis, 1981) and suggests that lava flow with a typical thickness (≫10 -1 m) takes >7 hr to solidify. Eruptive temperature is assumed as T e = 1,168 °C, and liquid viscosity is 64 Pa s (Soldati et al., 2021). The dashed black curve (1,100 °C) indicates the onset of viscosity increase to 218 Pa s (Soldati et al., 2021). The growing liquid-rich layer thickness is calculated by multiplying the time and estimated velocities from Equations 17-19. If the bubble size is relatively large (>0.01 m), bubble separation easily occurs, but this is not the case when the bubble size is small <1 mm. It is reasonable to consider that the internal bubbles have already experienced the coalescence and have a large size at the end of a lava flow.
It is helpful to consider a dimensionless number indicating whether bubble separation or solidification by cooling is more dominant. The thickness ratio of the liquid-rich layer L l and the cooled upper crust L c makes a dimensionless number. The evaluation of liquid-rich layer growth depends on the specific situation (Equations 17-19), but all cases depend on v stk , as defined by Equation 17. We here use v stk and neglect some constants for simplicity, and thus obtain
We note l∕ c ∝ √ , and bubble radius R will increase with time. Thus, Equation 21indicates that, if the flow thickness is sufficiently thick not to solidify too fast, the effect of bubble separation becomes more dominant with time. In a real lava flow, bubbles accumulate beneath the solidified upper crust which prevent growth of cold crust. We thus conclude that bubble separation can occur to form a cavity in a shelly paho'eho'e.
In Figures 4c and 5a, we observed the resuming of the stopped flow by the vertical separation of bubbles. This phenomenon is caused by the balance between the gravitational driving force of flow ρ l (1 -ϕ b )gh sin θ and the yield strength τ y . For highly crystalline magma, the yield strength will be high so that flow is likely to stop (Castruccio et al., 2014). In this case, the effective viscosity of the magma is also high, and bubble separation takes a long time. High bubble fraction also causes yield strength, but its value can be lower than that induced only by crystallinity (Birnbaum et al., 2021;Hohler & Cohen-Addad, 2005). Field measurements of yield strength of basaltic lava flow are in the order of 10 2 -10 3 Pa (Chevrel et al., 2019). Assuming a yield strength of τ y ∼ 10 3 Pa, melt density of ρ = 2,700 kg m -3 , vesicularity of 0.7, and slope angle of 25°, a 0.3 m thick flow would stop. Figure 11 suggests bubble separation can occur before the solidification in lava flows with a thickness of 0.3 m. We conclude that stopping and resuming of flow by changing the microstructure-dependent yield stress can occur in a natural lava flow.
We conducted a series of experiments of analog lava flows using bubbly corn syrup suspensions at a regime dominated by bubble surface tension (the capillary number, Ca < 1), where bubbles are expected to increase the effective viscosity of the suspension compared to the liquid viscosity. This regime is relevant for slow moving natural lava flows such as the flow end of thin overspills from main channels, and can work to slow the flow further than cooling alone. We identified a possibly non-Newtonian behavior during the early stages of the flow of bubbly suspensions, but then a return to Newtonian behavior. Liquid segregation and flow localization occurred in sections where bubble-poor, and thus a denser and also less viscous (since Ca < 1) fluid advanced faster than the rest of the flow. In natural lava flows, Ca > 1 can occur, which can accelerate flow velocity and enhance bubble coalescence. Large bubbles and liquid-rich regions can accumulate at the flow front. A bubble-rich region generated by the vertical separation might be an origin of shelly paho'eho'e.