The Yih Instability in Layered Lava Flow May Initiate the Pāhoehoe to ‘a‘ā Lava Transition

Most basaltic lavas begin flowing as pāhoehoe but sometimes transition into ‘a‘ā. Field observations and previous models have demonstrated that the rheology of smooth, liquid‐like pāhoehoe is distinct from rough, pasty ‘a‘ā, but the cause of this dramatic and rapidly occurring change in rheology has remained unclear. Here, we propose that the pāhoehoe to ‘a‘ā transition could be initiated by the Yih instability, an internal shear instability, in layered pāhoehoe flow. We show that the conditions under which instability arises depend on both extrinsic and intrinsic factors and derive a non‐dimensional regime diagram. We test the model's prediction of stable flow configurations against field observations of solidified lava flows.

Here, we propose that a mechanical stirring event can be triggered by the onset of the Yih instability, an internal shear instability that takes place in layered pāhoehoe flow. Our work is motivated by the finding that remolten pāhoehoe and 'a'ā both yield 'a'ā when stirred in the laboratory (Emerson, 1926), suggesting that stirring could play a role in the transition. Our contribution to this debate is to identify and characterize the specific physical process, an internal shear instability, that could lead to a sudden mechanical stirring and associated crystallization (Kouchi et al., 1986). By attempting to identify a specific physical process that could trigger the onset of the transition, our study complements field-based approaches to characterize the transition (e.g., Cashman et al., 1999;Macdonald, 1953;Peterson & Tilling, 1980;Soule et al., 2004;Swanson, 1973). To test our idea that layered pāhoehoe flow might be prone to dynamic instability, we devise a simplified mechanical model of a stratified pāhoehoe flow.
Our proposed model set up is shown in Figure 1. Shortly after exposure to the atmosphere, lava has approximately homogeneous flow properties, but exposure to the cool atmosphere soon forms a congealed upper surface ( Figure 1a). We consider the limit where the congealed-surface flows much slower than the layered lava underneath ( Figure 1b). To understand the stability of the flow underneath the congealed crust, we use a linear stability analysis (Yih, 1967) over a wide parameter range that is appropriate for pāhoehoe lava flows that are close to the "melting range" (Shaw, 1969). The linear stability analysis tests whether the presence of inevitable small perturbations η, in an otherwise stable density stratification, can exponentially grow with time t at a given location along a flow.
Our analysis draws on prior work by Yih (1967) and Yiantsios and Higgins (1988), who showed that stably layered viscous fluids can become unstable, forming an internal shear instability in space. We refer to this instability as the Yih instability to highlight that it is different from other, inertia-driven, shear instabilities such as the Kelvin-Helmholtz or Holmboe instabilities that commonly occur in atmospheric layers, but not in viscous lava flows. For example, the Kelvin-Helmholtz instability is visible in wave formations along clouds (Fritts & Rastogi, 1985), the Great Red Spot on Jupiter's atmosphere (Marcus, 1988) and many other wave-like features along low viscosity fluids (e.g., Burkholder et al., 2020;Woods, 1968).
A key assumption in our model is that pāhoehoe flows may have internal layers. Field observations of frozen pāhoehoe cross-sections at Deccan Traps (Aubele et al., 1988) and Western New Mexico (Duraiswami et al., 2014) support the possibility of multiple internal layers between the hotter inner lava and the cooler outer crust. Additionally, Cashman et al. (1999) argues from lava breakout temperatures that thermal stratification should be a critical component of any lava flow model. The existence of layers would also be compatible with the inflation of pāhoehoe sheet flows over areas that can be several hundred meters across (Hon et al., 1994) and domed structures called tumuli (Walker, 1991). Further evidence in support of layers are the small undulating surface structures along pāhoehoe lava flows that have been interpreted as the result of internal layering (Fink & Fletcher, 1978). Even in lava lakes, such as Lake Kivu at Nyiragongo volcano, Democratic Republic of Congo, surprisingly sharp differences in temperature have been documented at different depths within the lake (Lorke et al., 2004), despite the presence of convection which would tend to eliminate layering.
In addition to this observational evidence, layering in lava flows is also expected from a theoretical point of view. From emplacement to solidification, lava flows span a wide range of temperatures, volatile contents and crystallinities. The flow properties of lava, like effective density and viscosity, depend sensitively on these parameters. For example, temperature measurements from an active lava channel at Kılauea, Hawai'i, suggest temperature variations of about 200°C within the flow. This temperature difference could translate into viscosity differences of up to two orders of magnitude for melt with low water content, according to existing viscosity measurements (Hess & Dingwell, 1996). The loss of volatiles during conductive cooling and the presence of bubbles or crystals would tend to further exacerbate this effect (Costa et al., 2009;Lipman et al., 1985;McBirney & Murase, 1984;Shaw, 1969).
Prior research mostly models lava flows as a single layer with effective properties that depend on temperature, crystallinity, and bubble volume fraction (Cappello et al., 2016;Hidaka et al., 2005;Soule & Cashman, 2005). More recent models are now incorporating a congealed crust to enhance lava flow predictions (Cashman et al., 2006;Kilburn, 2004;Vicari et al., 2007), but retain the assumption of a single, relatively homogeneous layer of lava underneath. For the purpose of developing a general model for lava flow, considering the possibility of layering may constitute an unnecessarily complication, but the model we develop here is not intended as a general or realistic description of pāhoehoe flow. Instead, we are interested in understanding whether there are physical processes associated with layering that can help us understand specific aspects of lava flows such as the pāhoehoe to 'a'ā transition.
Here, we consider a layered, Poiseuille, pāhoehoe flow, where the top layer is less dense and more viscous than the bottom layer (see Figure 1b). We assume that both layers have a Newtonian rheology, implying that small bubbles and crystals are mixed in with the melt (Anderson & Jackson, 1967;Bowen, 1976;Slattery, 1967). This assumption limits our model to scales much larger than the bubbles and crystals, which is why we focus only on long wavelengths in our analysis. Following Yih (1967), we find steady, unidirectional solutions to the Navier-Stokes equations for each layer with appropriate boundary and interfacial conditions (see our Supporting Information S1, Equation 42 in Yih (1967) and Equation 11 in Yiantsios and Higgins (1988)). We then introduce an infinitesimal, long-wavelength perturbation at the interface between the two layers and analyze whether it grows or decays using linear stability as derived in Yih (1967) and detailed in our Supporting Information S1. Decaying perturbations imply that the flow remains stable unless it is subject to a major external disruption, while growing perturbations suggest the possibility of sudden flow instability. We set a no-slip boundary condition to model a surface speed that is much slower than the internal speeds and we assume a sharp boundary condition, such that the interface thickness is much smaller than the thickness of the flow.
Several non-dimensional terms arise in our linear stability analysis including the viscosity ratio m, layer thickness ratio n, density ratio r, and the extended Froude number Fr e . The extended Froude number represents the ratio between the kinetic and potential energy of the system. It is commonly used to describe geophysical mass flows on inclined slopes like avalanches and debris flows (Takahashi, 2007). Here, we define the extended Froude number squared as the ratio of the kinetic to the potential energy density in the upper flow layer: where d u is the upper-layer thickness and Δρ = ρ l − ρ u is the density difference between the layers. The characteristic speed, U o , could be related to the emplacement mechanism itself, gravitational effects or pressure gradients particularly in confined lava flows where a congealed-surface crust reduces the ability of the flow to dynamically adjust thickness.
We emphasize that the extended Froude differs conceptually from the way the Froude number is often defined in hydraulics to analyze hydraulic jumps. The Froude number in hydraulics has a characteristic speed that is governed by inertial forces, but relative to viscous forces, inertia is negligible for the flows considered here. Although inertial forces are negligible, the kinetic energy carried by the flow is not. This conceptual difference is reflected in the densimetric component of the extended Froude number. When ρ u ∼ ρ l , little potential energy is stored in the flow and kinetic energy always dominates, resulting in a high extended Froude number independent of flow speed.
We assume a viscosity difference of m = [0.1-0.33] between the layers, reflecting a <10 vol% crystallinity increase due to cooling in the upper layers . We impose a density ratio of 10% or r = [0.9-1], which could be the result of ∼10 vol % bubble increase in the upper layer. We test flow speeds that range between 0.01 and 30 m/s (Baloga et al., 1995;Lipman & Banks, 1987;Robert et al., 2014) and pick an arbitrary layer thicknesses d u = 1 m that is within field observations of 0.1-6 m (Aubele et al., 1988;Duraiswami et al., 2014) for the purposes of creating the regime diagram in Figure 2. The resulting ranges for the non-dimensional numbers we evaluate are n = [0.2-5] and Fr = [2.2 × 10 −2 −71].
In agreement with experimental studies (Miesen et al., 1992), the initial growth rate of instability scales with the wavenumber and the Reynolds number, which is the ratio of inertial forces to viscous forces, where the resulting range is . That said, given our mixture approximation, we do not consider small wavelength limit, which would be necessary to identify a fastest growing wavelength (Yiantsios & Higgins, 1988). Our analysis shows that pāhoehoe flow can become unstable for all large-wavelength perturbations if certain conditions are met. We quantify how relative viscosity, layer thickness, density, and flow speed between the layers impact the parameter space of instability in Figure 2. The regions between the dark blue and dark red, colored by light and dark salmon colors, are parameter regions where the stability of the interfacial perturbation is dependent on the viscosity ratio m between the two layers. Evidently, even relatively minor changes in the viscosity ratio, at least by volcanic standards, can alter the onset of instability notably.
We conclude that both extrinsic and intrinsic parameters affect instability. For example, flows with a large extended Froude number (i.e., flows with small density difference like Δ > 10 , high speeds greater than U 0 > 10 m/s, or thin upper layer thickness d u < 5 m), a thick lower layer (n > 1) and a high viscosity ratio (m < 1) are most prone to instability. Also, increasing flow speed U 0 , either through an increase in lava discharge or increasing slope inclination of the flow, will increase the parameter space of instability (Figure 2), because this configuration increases shear at the interface. We emphasize that the instability results summarized in Figure 2 could occur at both low or intermediate Reynolds number. Additionally, instability can still occur when the upper layer is more buoyant than the lower layer, which acts to stabilize instability.  Figure 2 shows a potential parameter space where the instability is applicable for the pāhoehoe to 'a'ā transition. However, demonstrating that a flow can become unstable is not the same as arguing that it will, which is important for determining if this instability does take place in lava flows. We examine field observations to determine if there is evidence to refute the applicability of our model for this instability. Parameter spaces that give rise to instability would not be preserved in layered field outcrops. Therefore, we can test the model prediction of stable (blue region of Figure 2) flow configurations against solidified pāhoehoe flows in the field, for example, the pāhoehoe flows that have cooled, and effectively preserved some flow characteristics like relative layer heights. We note that this comparison is inevitably flawed for a number of reasons, starting from the fact that we only consider the onset of instability but not model the dynamic evolution of the flow in response to it. Also, we only model pāhoehoe flows in the Newtonian limit, which limits our analysis to relatively small bubble and crystal fractions, while solidified pāhoehoe flows in the field have obviously gone through the entire evolution from largely fluid to solid. Nonetheless, we believe that it is essential to test model predictions against field observations, even if that comparison is imperfect.

Comparing Model Predictions Against Field Observations
We test our model against data provided in Aubele et al. (1988) and Duraiswami et al. (2014). Aubele et al. (1988) measured the relative layer heights of 29 solidified lava flows and flow fields in several volcanic fields erupted within the past million years in western New Mexico. We also integrate more recent data by Duraiswami et al. (2014), who used similar methods as Aubele et al. (1988) to measure an additional seven flows emplaced in the Bushe and Poladpur Formations of the Deccan Traps in India. They document three distinct zones in the vertical sections of these lava flows that are <16 m thick: an upper vesicular zone, a middle nonvesicular zone and a lower vesicular zone. This frozen flow stratigraphy with heavy vesiculation in the top and bottom layers suggest that each layer was part of a single lava flow that was exposed to cold temperatures along the top and bottom boundaries. Although we cannot be certain, this suggests that the layers were not multiple separate lava flows or thermally reactivated lava flows (Katona et al., 2021).
For our comparison, we neglect the presence of the lower vesicular zone, such that the middle nonvesicular zone thickness is the thickness of the lower layer in our model. A three-layered flow version of the Yih instability shows that the stability region stays comparable to the formulation we show here when the middle layer is less viscous than the top or bottom layer (Than et al., 1987). Although the layer thicknesses in the solidified flows can be observed directly (e.g., Aubele et al., 1988;Duraiswami et al., 2014), the density, viscosity, and speed of the two layers during active flow conditions are inevitably unknown. To compare the modeled and observed layer thicknesses, we sample density (r), viscosity (m), thickness (n) ratios, and lava flows speeds (U 0 ). We choose these limits based on observational constraints and model applicability. Aubele et al. (1988) reported that the bubble volume fraction in the upper vesicular zone is a few percent higher than in the nonvesicular zone; the dimensional flow speed of pāhoehoe tends to vary over a wide range of at least three orders of magnitude, for example, 0.01-30 m/s (Baloga et al., 1995;Lipman & Banks, 1987;Robert et al., 2014); and layer thicknesses may span 0.1-6 m, according to Aubele et al. (1988) and Duraiswami et al. (2014). These uncertainties lead to a wide span of possible Froude numbers characterizing the field sites. We adopt a probabilistic test to quantify the likelihood of the field data being consistent or inconsistent with the model. We choose to sample 10,000 random sets of parameters to introduce variability in our sampling. Any sampling greater than 10,000 changes the median layer thickness ratio results by less than 10%. There is significant overlap in the stable and unstable layer thicknesses, because the layer thickness ratio is not the only determining factor of this instability. Our model's median height ratio for stable flows, n = 0.84, lies similar to the medians of the pāhoehoe field studies, n = 0.64 (Duraiswami et al., 2014) and n = 0.84 (Aubele et al., 1988). Meanwhile, the median ratio for unstable flows is n = 1.12, which is still within the standard deviation of the stable flow height ratios. Additionally, Figure 3 shows the number N of stable and unstable flow configurations, which add up to 10,000, the number of random initializations. Our results show that most modeled flows are stable and only 10% of the cases for the given parameter range become unstable.
The comparison shown in Figure 3 is by no means a validation of our model and not intended that way. We perform this comparison to field data to refute our model, but do not find any major inconsistencies between our model results and the field observations of Duraiswami et al. (2014) and Aubele et al. (1988). We hence do not see grounds to reject the model. Narrowing the scope of any one of the parameter ranges can notably decrease the overlap, but there is no strong justification for narrowing these parameters for the given field sites. We cannot conclusively rule out the possibility that the data fall into the range of unstable layer thicknesses.

Discussion
There are multiple processes that can contribute to the pāhoehoe to 'a'ā transition. Prior research has proposed the transition could be due to mixing following lava falls, but lava downstream of lava falls may still flow as pāhoehoe (Cashman et al., 1999). The transition could also be due to outbursts of stagnant lava that has been cooling underneath a congealed crust (Cashman et al., 1999;Peterson & Tilling, 1980). This idea is consistent with the observation that the transition can appear to be sudden (Peterson & Tilling, 1980), but is inconsistent with the observation that the transition has a correlation with high flow rates (Macdonald, 1953;Peterson & Tilling, 1980;Pinkerton & Sparks, 1976;Rowland & Walker, 1987, 1990. Finally, the transition could be due to the high strain rates near the spontaneous formation of "clots" along pāhoehoe flows (MacDonald, 1958;Peterson & Tilling, 1980;Wentworth & Macdonald, 1953). These clots possibly form due to portions of the flow surface being introduced into the flow interior (Kilburn, 1990(Kilburn, , 2004.
None of these ideas are incompatible with the possibility we are considering here, namely that the Yih instability may trigger the initiation of the pāhoehoe to 'a'ā transition. External disruptions like lava falls undoubtedly perturb the internal flow field, but whether this perturbation decays or grows depends on whether the internal flow field is close to instability or not. An outburst of lava from stagnant flows results in an increase in internal shearing, which could trigger the transition when the intrinsic conditions are conducive to instability. Alternatively, tearing of the congealed surface and the sudden introduction of lava clots into the flow constitute the types of inevitable disruptions that could trigger the onset of instability when the flow is on the cusp of instability. Our results are hence complementary to these existing ideas, adding a potential explanation for why external factors alone do not appear to fully control the onset of instability. If the flow is sufficiently far from becoming unstable, even a significant disruption might not be enough to initiate instability, but when the flow is on the verge of becoming unstable, even a very subtle change in conditions could initiate instability.
While the dependence of the pāhoehoe to 'a'ā transition on both intrinsic and extrinsic properties has certainly been pointed out before (Soule et al., 2004), our work quantifies this dependence through the regime diagram shown in Figure 1. Intrinsic properties characterize the lava itself, such as composition and volatile content, and are represented in our model through their effect on lava density and viscosity. In contrast, extrinsic properties characterize the flow, such as its slope and distance from the vent, and are represented in our model through lava flow speed (Figure 1b). The regime diagram in Figure 2 suggests that a wider range of intrinsic properties are associated with instability for fast as compared to slow flowing pāhoehoe (e.g., high Fr). The finding that fast flow is more prone to instability is consistent with the observation that the transition is often (Cashman et al., 1999;Macdonald, 1953) but not always Soule et al., 2004) associated with an increase in the slope of the underlying terrain. If the flow is sufficiently far from becoming unstable, even a significant change in one parameter such as slope might not be enough to initiate instability, explaining why a change in slope matters sometimes but not always. Similarly, a very subtle change in conditions could initiate instability if the flow is already on the verge of becoming unstable.
The possibility that the Yih instability contributes to the initiation of the pāhoehoe to 'a'ā transition may explain why the pāhoehoe to 'a'ā transition appears to be mostly, but not strictly, irreversible (Kilburn, 2000;Peterson & Tilling, 1980). It is common for flows to transition from pāhoehoe to 'a'ā, but the inverse transition from 'a'ā to pāhoehoe is rare Kauahikaua et al., 2003). The mechanical shear instability we focus on here is irreversible, because it is associated with an increase in entropy. Yet, recent observational literature conclude that the process could be reversible in rare cases, especially for long-lived eruptions of low and/or variable effusion rates Kauahikaua et al., 2003). Within our model framework, the reverse transition from 'a'ā back to pāhoehoe lava could be a consequence of the introduction of new, hot, low viscosity lava into the lower layer in Figure 1b. The hot lava then either breaks out through the 'a'ā, giving the appearance of the flow transitioning back, or remelts the crystals in the 'a'ā, facilitating an actual transition back.
An important caveat of our model is that it assumes a congealed surface has formed. Shortly after eruption, lava flows resemble free-surface flows with approximately homogeneous properties (see Figure 1a). Perturbations along a free surface are also subject to instability, which has been studied analytically (Benjamin, 1957;Yih, 1954Yih, , 1963 and experimentally (Binnie, 1957;Kapitza, 1948;Liu et al., 1993). Experimental research and analytical work show that flow angle and Reynolds number govern the onset of the free surface instability, where there is a bifurcation point at Re = 5/6 cot(θ) that characterizes the internal (Re) and external (θ, the inclination angle) conditions when an instability occurs (Liu et al., 1993;Yih, 1963). These results suggest that steep and fast lava flows characterized by Re greater than 0.01 could be subject to this free-surface instability.
Although free-surface instabilities could be relevant in the immediate vicinity of an eruption vent, rapid crust formation due to cooling hinders their growth rate and limits their practical importance for active lava flows. Maybe more importantly, observations show that the pāhoehoe to 'a'ā transition does not always tend to occur in the immediate vicinity of vents (Duraiswami et al., 2014;Peterson & Tilling, 1980;Soule et al., 2004). Thus, a potential important dynamic role of the crust is to incubate the lava flow underneath and thereby enable a greater time period for the Yih instability, highlighting the value in investigating the flow behavior between the limits of a free-surface and congealed-surface flow. Given the subtle interplay between free surface and internal dynamics and the possibility of multiple instabilities, a numerical solution to the instability problem would then become necessary.
In addition to considering surface-driven perturbations, lava flows may be influenced sensitively by mesoscale topographic wavelengths in the range of 10-100s of meters (Richardson & Karlstrom, 2019). The model we propose here could be extended to better understand the impact of basal topography by adding a basal layer to our model setup. As suggested by the observations of Aubele et al. (1988), this lower layer would likely be cooler and more vesicular than the layer above and might act to filter out some, but not all of the perturbations inevitably created by flow over basal topography. Again, a generalized numerical or experimental investigation would be necessary to characterize the implications of basal topography as a perturbation source for the Yih instability. In the interest of providing a first step toward a more complete understanding, we intentionally adopt a simplified model that can be solved analytically and advance our process-based understanding for why abrupt changes in flow behavior may occur.
Our model only identifies the onset of instability, but linear instability does not necessarily imply nonlinear instability. Therefore, our analysis merely identifies the onset of instability but does not apply once the instability has reached a finite amplitude. Even in two completely homogeneous layers, nonlinear processes like viscous damping of finite-sized waves and thermal diffusion will alter the progression of instability soon after onset.
In multiphase layers, crystallization and bubble formation may enhance or hinder instability growth. Given that instabilities take time to grow and there is inevitable heterogeneity in the evolving mixing process, we would expect the transition to have a clear onset, but our model does not constrain what the post transition flow or rheology may look like. The post transition flow could become metastable where nonlinear effects suppress continued growth of the instability, often resulting in standing waves on the interface, or the flow could become unstable, likely leading to a sudden mixing of the two, initially separate layers. The post transition rheology may resemble a more viscous pāhoehoe, like those described by many observational studies (Duraiswami et al., 2014;Jones, 1943;Malin, 1980;Murcia et al., 2014;Rowland & Walker, 1987) or may include both flow types coexisting, like that described in Soule et al. (2004) before the transition completes.
One potential example of this instability where the transition may look different is at the Nyiragongo volcano, DRC. Nyiragongo produces ultra mafic lava flows (Aoki et al., 1985) with initial lava viscosities that can be as low as 60 Pa·s (Giordano et al., 2007), which should travel rapidly over long distances. However, Nyiragongo has steep edifice, which is characteristic of more viscous lava. Lava at Nyiragongo flows rapidly until it experiences rapid crystallization and slows down, resulting in the formation of steep edifices (Morrison et al., 2020). These lava flows may be experiencing rapid crystallization because of the Yih instability. Our finding that fast flowing lava over steep terrain is particularly prone to the Yih instability could explain the rapid crystallization and formation of a steep edifice ( Figure 2).
Here, rather than modeling the pāhoehoe to 'a'ā transition, we merely propose a specific physical process that could control its initiation. Identifying a specific physical process, in our case the onset of the Yih instability, is different from identifying a single physical parameter like surface slope. The shift we propose, from a specific governing parameter to a specific process, is motivated in part by prior research that has struggled to clearly identify a single parameter governing the transition (Murcia et al., 2014;Peterson & Tilling, 1980;Soule et al., 2004). We make an attempt to test our model results against solidified lava flows while keeping in mind that this comparison is inevitably flawed. The next step would be to compare the model to active lava flows as documented in Dietterich et al. (2022) and Cashman et al. (1999), ideally once the model has been generalized to include temperature and nonlinear effects. With the deployment of new technologies that track activity in Kılauea, Hawai'i and other volcanoes (Thompson & Ramsey, 2020) and more direct observations of active lava flows (e.g., Dietterich et al., 2022), new opportunities might soon arise to conduct further, more rigorous testing of our ideas.

Data Availability Statement
We provide the scripts used to solve the matrices for the system of equations and create the plots at our Zenodo repository here: https://doi.org/10.5281/zenodo.7478513. The data used for Figure 3 can be found in Aubele et al. (1988) and Duraiswami et al. (2014).