Research ArticleGEOPHYSICS

Crystal aggregates record the pre-eruptive flow field in the volcanic conduit at Kīlauea, Hawaii

See allHide authors and affiliations

Science Advances  04 Dec 2020:
Vol. 6, no. 49, eabd4850
DOI: 10.1126/sciadv.abd4850


Developing reliable, quantitative conduit models that capture the physical processes governing eruptions is hindered by our inability to observe conduit flow directly. The closest we get to direct evidence is testimony imprinted on individual crystals or bubbles in the conduit and preserved by quenching during the eruption. For example, small crystal aggregates in products of the 1959 eruption of Kīlauea Iki, Hawaii contain overgrown olivines separated by large, hydrodynamically unfavorable angles. The common occurrence of these aggregates calls for a flow mechanism that creates this crystal misorientation. Here, we show that the observed aggregates are the result of exposure to a steady wave field in the conduit through a customized, process-based model at the scale of individual crystals. We use this model to infer quantitative attributes of the flow at the time of aggregate formation; notably, the formation of misoriented aggregates is only reproduced in bidirectional, not unidirectional, conduit flow.


Located on the Big Island of Hawai’i, Kīlauea volcano threatens hundreds of homes and millions of dollars of property. To detect potential unrest and eruptive activity, Kīlauea is constantly monitored by the Hawaiian Volcano Observatory through an extensive array of field stations, recording seismic, geodetic, acoustic, and gas-flux data. The resulting, rich datasets provide a unique opportunity for improving theoretical models of volcano dynamics, but it is not obvious how to test models of physical processes occurring at depth against observations recorded at the surface. As a consequence, we still lack reliable quantitative models describing the flow regime in the volcanic conduit, which limits our ability to predict eruptive activity and mitigate the risks it poses.

Geophysical observations provide constraints on the physical processes operating at the scale of the volcanic edifice and the surrounding crust. Examples include deformation measurements of the volcanic edifice over time (1, 2) and seismic imaging of its subsurface structure (3, 4). While invaluable for understanding the large-scale dynamics, geophysical observations are less suitable for identifying the flow regime inside the conduit. In contrast, crystals trapped in magma quenched upon eruption record at least some aspects of the pre-eruptive condition in the volcanic conduit directly. Currently, crystal-scale data are most often used to estimate the ascent rate (5) that is then used as an input to models. However, we rarely attempt to test the underlying physical assumptions of different existing conduit models against crystal-scale data. The main reason is that most conduit models operate at the scale of meters or kilometers (6, 7) and hence do not entail testable model predictions at the scale of individual crystals. As a consequence, we may not fully leverage the potential treasure trove of information imprinted onto individual crystals in the volcanic conduit.

Here, we link two crystal-scale models to extract constraints on the pre-eruptive conduit flow from small olivine aggregates erupted in picritic scoria at Kīlauea Iki in 1959. Examples of these olivine aggregates, analyzed by Schwindinger and Anderson (8), are shown in Fig. 1 (A to C). Within a given cluster, olivine crystals orient preferentially along certain crystallographic faces (see Fig. 1D). Most notable is the substantial percentage (37%) of {021}-{021} aggregates that feature a large angle between crystals (Fig. 1E), commonly referred to as a misorientation angle (9). Large misorientation angles are unexpected because, relative to two parallel crystals, they imply a large surface area that elevates hydrodynamic drag. The common occurrence of aggregates with large misorientation angles hence requires a mechanism that creates this particular orientation (810). The physical origin of this mechanism, however, has remained elusive. Previous attempts to recreate the observed orientations through crystal settling in the laboratory have not been successful (11, 12).

Fig. 1 Crystal observations, map, and photo from 1959 eruption at Kīlauea Iki.

(A to C) Selected examples of olivine aggregates found in Kīlauea Iki scoria, Iki 22, S-5 (27) from Fig. 2 in Schwindinger and Anderson (10). The crystals contain inclusions of chromite (black), glass (gray), and gas [black circles with white spots, as in the upper left corner of (B)]. (D) Crystallographic faces of the olivine crystals from (10). (E) Summary of alignment proportions as measured by (10), compiled by K. Allison. (F) Map of the Kīlauea summit area and lava flows generated by the 1959 eruption, redrawn after Fig. 3 in Richter et al. (26) by K. Allison. (G) Location of Kīlauea volcano on the Big Island of Hawaii. (H) Photo of a lava fountain blasting lava fragments 425 m high above the vent during the third fire-fountaining episode of the eruption. Photo credit: U.S. Geological Survey.

We hypothesize that the large misorientation angles in the Kīlauea Iki aggregates were produced by exposure of the olivine crystals to a steady wave field in the conduit. Our hypothesis is motivated by the documented ability of wavy flows to align particles (13, 14) in a manner distinct from the behavior of crystals in shear flow. Of the two most common classes of conduit models, positing unidirectional [see Sahagian (15) for a recent review and a comprehensive list of references] and bidirectional flow (1623), only bidirectional models naturally entail a steady, traveling wave (see Fig. 2). The wave originates at the interface between buoyant, gas-rich magma ascending in the center of the conduit and denser, degassed magma descending on the sides (24). It is a reflection of the linear instability of the interface (25). The growth of the instability is suppressed by a finite viscosity contrast between ascending and descending magma, leading to a metastable flow configuration with a stable interface wave (18, 22, 24).

Fig. 2 Schematic of the two classes of conduit models considered.

Illustrations, not-to-scale, represent unidirectional (left) and bidirectional (right) flow and the flow field they entail at the crystalline scale. The model domains used in our modeling are highlighted in light gray. The model domain for unidirectional flow is linear shear and for bidirectional flow it implies an interfacial wave superimposed on a linear shear flow at the interface, x = 0. The latter model domain is shown in the reference frame traveling at U0, so there is no mean flow at the interface. In both cases, the shear is characterized by the gradient du/dx¯=γ. u, velocity; x, horizontal direction; z, vertical direction.

To test our hypothesis, we quantify crystal alignment and aggregation in a pure-shear (unidirectional) flow, a pure-wave flow, and shear flow imposed by a steady, traveling interfacial wave (bidirectional flow) (see Fig. 2). A pure-wave field is unlikely in volcanic conduits, but we include it here for comparison purposes. We first derive an analytical criterion to identify the condition under which crystals align under the joint influence of wave and shear flow, building on DiBenedetto et al. (13). Second, we conduct direct numerical simulations based on the methodology by Qin et al. (12) to quantify the relative orientation of multiple interacting and eventually colliding crystals with respect to each other. Throughout the paper, we use the term alignment to describe the angular relationship between a single crystal and the main flow direction in the ambient flow field. In contrast, we use the term orientation to describe the angular relationship between two or multiple crystals with respect to each other (see Fig. 4, C to D) (9).

Together, our analytical and numerical models demonstrate that small crystal aggregates with large misorientation angles as observed in the Kīlauea-Iki scoria samples from 1959 form naturally in bidirectional but not in the more commonly used framework of unidirectional conduit flow. In addition to distinguishing the degree to which these two conduit models are compatible with crystal-scale observations, our models allow us to draw several quantitative inferences about the conduit flow field at the time of aggregate formation, including the relative magnitudes of wave- and shear-forcing, the direction of the wave, the spatial distribution of crystallinity in the conduit, and the likely location in the conduit where the samples originated. These constraints could shed light on the physical processes leading to eruption and why the samples were erupted as scoria rather than liquid lava during fire-fountaining.

The Kīlauea Iki Glomerocrysts

The 1959 eruption started on November 14 in Kīlauea Iki, Hawaii, a collapsed crater adjacent to the main summit caldera at Kīlauea (see Fig. 1, F and G). Overall, the entire summit eruption entailed 17 different fire-fountaining episodes and ended in March 1960 after the collapse of the Halema’uma’u crater (26). Figure 1H shows a photo of a lava fountain during the third eruptive episode. Apart from the incandescent magma feeding the lava flows, a small portion of magma is erupted as pumice or scoria (black) and blown off to the southwest, accumulating into the pumice blanket and forming scoria spatter cones (see Fig. 1F).

The relative accessibility of the eruption allowed systematic sampling of the outpouring magma to document chronological changes in composition (27). After pairs of samples collected from both the lava flows and the scoria had indicated no notable differences in silica content (27), only scoria samples were collected twice daily during the remainder of the eruption. Here, we consider the same scoria samples as Schwindinger and Anderson (10), erupted mostly during the first fire-fountaining episode of the summit eruption with one sample each from the third and eighth episodes.

Most aggregates contain between two and four crystals, but they can be as large as 16 crystals on rare occasions (see Fig. 1, A to C) (10). About half (54%) of all olivine crystals are part of an aggregate (11). Within an aggregate, most individual crystals are roughly comparable in size and have aspect ratios around two (10). Figure 1E summarizes the proportion of crystals aligned along particular crystallographic faces for a total of 244 aggregates (10). The two most common orientations observed are small misorientation angles of approximately 0, characterizing both the {010}-{010} and {110}-{100} aggregates, which together constitute about 40% of all aggregates, and large misorientation angles of approximately 80 (9) observed in the {021}-{021} aggregates, which constitute about 37% of all aggregates. We note that these percentages differ from the later analysis by Wieser et al. (9), which included samples from several historic eruptions in addition to the scoria samples from 1959.

The Kīlauea Iki aggregates are thought to have formed through synneusis (9, 10), defined as the drifting together and attachment of crystals through long-range hydrodynamic interactions (28, 29). Several pieces of evidence support this interpretation including the similar size of aggregated crystals (10), the broad range of misorientation angles (9, 10), and electron microprobe analysis (8), indicating that some of the crystals appear to have originated in magma with different silica content and temperature (10). Variable proportions of glass inclusions, chromite, and zonation patterns also suggest that many of the crystals have different provenance, while some appear to have grown under comparable conditions and, potentially, close to each other (8, 10). Most of the contact zones between crystals are partially or completely overgrown, suggesting that the crystals were in contact for an extended period of time (10). Some contact zones have glass wedged in between the crystals (10).


Crystals align at finite and opposing angles with respect to the direction of a steady, traveling wave

Individual crystals behave differently in shear-dominated flow as compared to a wave-dominated flow (30, 31). To illustrate the difference, we first examine the dynamics of individual crystals in each flow in isolation before superimposing the two. In pure shear, an ellipsoidal crystal will continually rotate, undergoing what is referred to as a Jeffery orbit (32), as sketched in Fig. 3A. These orbits can be broken by crystal inertia, fluid inertia, crystal-crystal interactions, or changes to the flow (33, 34). In a pure wave, small crystals with finite aspect ratio assume a preferential alignment, which is not parallel to the flow direction. They exhibit a small wobble around the preferential alignment angle as shown in Fig. 3B. The crystal behavior represents a stable limit cycle about a fixed point, ϕ¯* (13, 35), where ϕ denotes the particle’s alignment with respect to the ambient flow field (Fig. 3C), ϕ¯ is the wave-period–averaged orientation, and * identifies a fixed point.

Fig. 3 Individual crystal orientation behavior under a shear flow and under a propagating wave.

(A) Crystal tracing out a Jeffery orbit in pure, upwards-oriented shear. (B) Crystal wobbling around the preferential alignment angle induced by a pure wave propagating upwards. (C) We define the alignment angle, ɸ, as the angle between the crystal’s axis of symmetry and the x axis, which points horizontally across the conduit. (D) We define the misorientation angle, θ, as the smaller angle between two crystals’ axes of symmetry in an aggregate. (E) The four possible crystal configurations of individual, preferentially aligned crystals with axis-symmetry. The red and gray crystals are aligned at ±ϕ¯*, respectively. The top two configurations represent large misorientation angles and the bottom two show parallel crystals. (F and G) Crystal orientation over time for a single particle in pure shear (F) as compared to a wave field (G). The light gray curves represent the results from the analytical model, and the black curves denote data from the direct numerical simulations. Time is normalized by the Jeffery period TJ (F) and by the wave period Tw (G), respectively. u, velocity; x, horizontal direction; z, vertical direction.

Figure 3 (F and G) compares the temporal evolution of the alignment angle ϕ in pure-shear (Fig. 3F) and pure-wave flow (Fig. 3G) for the analytical model (gray) and the numerical model (black). The crystals have the same aspect ratio in both models, but the specific crystal geometry is different. We assume a smooth, ellipsoidal shape in the analytical model and an angular, rectangular shape in the numerical model. Another confounding factor in the comparison is that we assume a linear wave in the analytical model while the waves in our simulations are inevitably nonlinear, if only mildly so. Despite these differences, the overall, qualitative behavior of the crystals is similar in the two flows (see Fig. 3, F and G), suggesting that the aspect ratio is the primary factor governing the angle at which axisymmetric crystals align when keeping wave conditions constant. The specific number of crystallographic faces of the crystal and their respective angles is hence less important than the overall aspect ratio in determining alignment behavior.

To understand crystal orientation in bidirectional flow where both flow features, shear and waves, are superimposed and create competing effects on crystal alignment, we derive an analytical criterion for when single crystals exhibit preferential alignment. Our model is based on a Taylor expansion of the governing equations as detailed in Materials and Methods. We show that the transition between orbiting behavior and preferential alignment is a function of both the crystal aspect ratio λ and the relative strength of the mean shear rate γ compared to the shear rate associated with the wave transport γw, characterized by 𝒟 = γ/γw.

We find two clear limiting cases: If the flow gradient due to the waves is much stronger than the mean shear, ∣𝒟 ∣ ≪ 1, then individual crystals exhibit a stable limit cycle about a mean angle while undergoing small-amplitude oscillations at the wave frequency. In the limit where the mean shear is much stronger than the wave-induced flow gradients, ∣𝒟 ∣ ≫ 1, we see modified Jeffery orbits, where the tumbling orbit seen in the pure-shear case has a superimposed oscillation at the wave frequency. These findings demonstrate that the presence of a wave can substantially alter crystal behavior by stabilizing preferred crystal alignment, even in a flow field with shear superimposed.

In Materials and Methods, we derive the exact preferential alignment angle and find it to be a function of both the particle eccentricity ε, related to the aspect ratio as ε = (λ2 − 1)/(λ2 + 1), and the shear ratio 𝒟, whereϕ¯*=12cos1(ε+D/ε1+D)(1)

From this expression, it follows that the only real solutions for ϕ¯* exist for ∣𝒟/ε ∣ ≤ 1. For ∣𝒟/ε ∣ > 1, the preferential alignment breaks down and Jeffery orbits reemerge. As shown in Fig. 4A, as 𝒟/ε varies from −1 to 1, the alignment angle spans all possible values for one crystal. The time to approach the angle is not constant and depends on the initial condition. In Fig. 4B, the relationship between the preferential alignment angle ϕ*¯ and crystal eccentricity ε is shown as a contour plot. For each particle shape, from oblate (ε < 0) to prolate (ε > 0), the most extreme alignment angles (ϕ¯*=0,90) correspond to the strongest shear scenarios (𝒟/ε = − 1,1, respectively).

Fig. 4 Preferential alignment of a single crystal under wave and shear flow.

(A) Wave-period-averaged alignment ϕ¯ over time for a single crystal (ε = 0.6) in a mixed wave and shear flow with variable 𝒟/ε (see color bar on the right). Data come from numerically integrating the analytical governing equations (18, 20, and 19). The dashed line denotes 𝒟/ε = 0, or the case of only waves and no mean shear flow. (B) Contour plot showing predicted ϕ*¯ from Eq. 1. Color bar represents 𝒟/ε. The dashed line denotes the curve for 𝒟/ε = 0. The vertical yellow line corresponds to the flow conditions that would lead to misorientation angles observed for the Kīlauea Iki olivines. ε, crystal eccentricity; 𝒟, shear ratio.

The Kīlauea Iki olivines tend to have aspect ratios around two (λ = 2, ε = 0.6) and are hence prone to exhibit preferential alignment when exposed to a sufficiently strong, steady wave field. The axial symmetry of olivine causes the preferential alignment of the crystals to be a set of two angles, ±ϕ¯*, which correspond to a reflection about the x axis pointing in the horizontal direction across the conduit (see Fig. 3C). Thus, two crystals of the same shape can align onto two opposing angles in the same flow field. Because of these two alignment angles, four possible aggregate combinations exist as shown in Fig. 3E. If one crystal is aligned onto +ϕ¯* and the other is aligned onto ϕ¯*, then the crystals become misoriented when aggregated. The implied misorientation angles from these configurations are θ=2ϕ¯* or θ=1802ϕ¯*. If both crystals are aligned onto the same angle, then the crystals are parallel when aggregated.

These alignment combinations may provide a process-based explanation of the notable prevalence of the large misorientation angles in the 1959 Kīlauea Iki samples (see Fig. 1E). To identify the flow conditions consistent with observations, we use the misorientation angle of 80 estimated by Wieser et al. (9). This misorientation angle could result from a preferential alignment angle of either 40 or 50 (see Fig. 3E), highlighted as a vertical yellow bar in Fig. 4B. The marked region coincides approximately with 𝒟/ε ≈ − 0.5, indicating that a wave traveling in the opposite direction of the shear flow is necessary to reproduce the misorientation of the observed crystals with ε = 0.6. Since the flow direction in the ascending core is upward oriented, the wave we infer here must travel downward. Together, our analytical analysis hence suggests that the Kīlauea Iki olivines were aligned by a pronounced, steady, and downward traveling wave.

Previous work demonstrated that both upward and downward traveling waves can form on the bidirectional flow interface (22, 24). In the specific context of core-annular flow, two different flow configurations are equally viable in steady state, one associated with a thick, relatively slow-moving core and one representing a thin, much faster core flow (22). A downward traveling wave arises only when the radius of the ascending core flow is sufficiently large (22, 24) and may hence be indicative of the thick-core flow configuration more commonly observed in the laboratory (18, 21).

Finite crystal fraction suppresses large misorientation angles

Our analytical model only quantifies the effect of different flow fields on a crystal, but does not account for the reverse, namely, the effect that the crystal has on the flow field. As a consequence, it does not capture the long-range, hydrodynamic interactions between multiple crystals that arise as a consequence of the crystals’ joint effect on fluid flow. To test whether the wave-induced set of two preferential alignment angles remains robust when crystals are interacting and ultimately colliding, we use direct numerical simulations (12) to generalize our analytical findings to finite crystal fractions.

Figure 5 shows snapshots of three simulations at 5% crystal fraction in pure-shear flow (Fig. 5A), a mixed flow field with shear and a dominant wave component superimposed (Fig. 5B), and a pure-wave field (Fig. 5C). The mixed flow field is informed by the results of the analytical model, suggesting that a set of two alignment angles consistent with the observed misorientation angles emerges at 𝒟 ≈ − 0.5. This value represents a flow where the shear rate associated with the wave is approximately three times stronger than the mean linear shear rate, and the wave is traveling downward into the conduit.

Fig. 5 Direct numerical simulations of crystal motion over time for the three flow fields.

All simulations are initiated with a crystal fraction of 5% and crystals distributed randomly throughout the domain. The flow fields shown are pure shear (A), a mixed flow field with shear and wave superimposed at 𝒟/ε = −0.5 (B), and a pure wave field (C). Selected crystal contacts are marked in red. The colors associated with the three different flow fields corresponds to the colors used in Fig. 6.

A notable difference between the three cases shown in Fig. 5 is the heterogeneous distribution of crystals within the model domains. Crystals in pure shear remain almost homogeneously distributed throughout the domain (Fig. 5A). In contrast, crystals in the mixed flow field (Fig. 5C) and the pure-wave field (Fig. 5B) are displaced almost entirely into the half of the domain that is less exposed to the interfacial wave. These simulations demonstrate that the flow field affects both the spatial distribution of crystals and the alignment of individual crystals. The tendency of crystals to distribute heterogeneously is a consequence of the long-range, hydrodynamic interactions between individual crystals that, in the absence of inertial forces limiting their spatial range, break the symmetry of the flow in a similar way as turbulence in the inertial limit (36), as demonstrated in previous studies (12, 37, 38).

Figure 5 also shows that crystals in all three flow fields display a range of individual alignment angles with respect to the flow. In pure shear (Fig. 5A), a wide range of alignment angles is expected as individual crystals are undergoing Jeffery orbits in the flow. However, we also find a range of alignment angles in the two other cases (Fig. 5, B and C), where we would expect individual crystals to align preferentially, if not quite as wide a range as in pure shear (Fig. 5A). This finding demonstrates that the hydrodynamic interactions between crystals alter the wave-induced preferential alignment of individual crystals even at low crystal fractions. Contrary to the analytical model, where the preferential alignment angles depend on the propagation direction of the wave, we find no notable difference of the crystals’ alignment for waves propagating in different directions in our simulations. This difference is likely the consequence of the simulations being performed in two dimensions, where aspect ratio is the main descriptor of crystal geometry. We are hence not able to resolve differences in alignment due to variations in crystal eccentricity.

To better understand the implications of the wave-induced preferential alignment on the misorientation angles in crystal aggregates, we highlight crystals in contact through red circles in Fig. 5 (A to C). We emphasize that crystal contacts in our simulations are transient, because we assume that crystals do not stick upon initial contact. They only remain in continued contact if the hydrodynamic forces are amenable to that, otherwise they drift apart again. Similarly to individual crystals, contacting crystals exhibit a wide range of misorientation angles in the three different flow fields. At first sight, contacts with large misorientation angles appear to occur in all three cases. A closer look, however, reveals an important difference. In pure-shear flow, the individually rotating crystals collide with each other, but they keep rotating individually albeit at a slower rate even when in contact with another crystal (Fig. 5A). As a consequence, the colliding crystals remain at least partially mobile as apparent from the visible changes in misorientation angle for the two contacts highlighted at four different times in Fig. 5A. In contrast, contacts in a wave-dominated or pure-wave field are relatively stable, because crystals mostly remain in their original alignment after coming in contact (see Fig. 5, B and C).

At low average crystal fraction (5 to 15%), long-lived crystal contact occurs only in wave-dominated flow

The subtle differences in crystal contact and misorientation shown in Fig. 5 call for a more systematic, statistical analysis. To compare our simulations to field observations, we plot the distributions of misorientation angles (9) in Fig. 6. Since we simplify the shape of the olivines from polygonal prisms to ellipses in the analytical model and rectangles in the numerical model, we are not able to distinguish between individual crystallographic faces. We categorize both {110}-{110} and {010}-{010} contacts as parallel with a misorientation angle of less than 18. The relatively rare {021}-{010} contact corresponds to a misorientation angle between 36 and 54, and the {021}-{021} contact represents a large misorientation angle of more than 72. In Fig. 6, we show the observed percentage of aggregates with large misorientation angles (37%) (10) as a gray block behind the simulation results.

Fig. 6 Misorientation angles and longevity of crystal aggregates from direct numerical simulations.

(A to C) Histograms showing the distribution of misorientation angles during collisions in a pure-shear, a mixed, and a wave field at 5% crystal fraction (A), 10% crystal fraction (B), and 15% crystal fraction (C). Grey bars show 37% for comparison with observed samples. (D to F) Histograms showing the longevity of crystal contact after collision for all three wave fields at 5% crystal fraction (D), 10% crystal fraction (E), and 15% crystal fraction (F).

Crystal contacts in all three flow fields exhibit the full range of misorientation angles (see Fig. 6, A to C), highlighting the stochastic nature of crystal collisions in flow. The most common crystal configuration in all three flow fields is a parallel or almost-parallel contact as expected because this configuration minimizes surface area and hence minimizes hydrodynamic drag. The only case in which the percentage of contacts with a large misorientation angle is comparable to observations is a pure-wave or a wave-dominated flow field at 5% crystal fraction (see Fig. 6A). We note that we do not exactly reproduce the 37% of large misorientation clusters observed in Kīlauea Iki samples (8). Multiple factors may be contributing to the slight discrepancy including uncertainty in measurement, stochastic variability, classification of the observed misorientation angles, and the simplified crystal shape and waveform considered in our models. When averaging more than five simulations for each flow case, about 30 ± 5% of contacts in a pure-wave field and 26 ± 5% of contacts in a wave-dominated field with 𝒟/ε = − 0.5 have large misorientation angles. Even a slight increase in crystal fraction from an average of 5 to 10% suppresses the wave-induced alignment of crystals sufficiently to all but eliminate crystal contacts with high misorientation angles (see Fig. 6B).

One confounding factor when comparing the statistics of instantaneous contacts as shown in Fig. 6 (A to C) to the overgrown olivine aggregates found in Kīlauea-Iki scoria is the longevity of contact. In the field samples, most of the contact zones between the olivines are partially or completely overgrown, suggesting that the crystals were in contact for an extended period of time (10). In Fig. 6 (A to C), we plot the misorientation angles in instantaneous contacts, without considering the duration of contact. In Fig. 6 (D to F), we show histograms of the longevity of crystal aggregates for the three flow fields and crystal fractions. These data demonstrate that the longevity of contact also depends sensitively on the ambient flow field. In pure-shear flow, the vast majority (≈90%) of crystal contacts is extremely short-lived for the entire range of crystal fractions considered (see blue bars in Fig. 6, D to F). In a wave-dominated and pure-wave field, crystal aggregates are stable. For example, over half of all crystal contacts in wave-dominated or pure-wave flow last longer than 80 wave periods at an average crystal fraction of 5% (Fig. 6D). The longevity of contact declines slightly as the crystal fraction increases to 10 or 15% but remains elevated at about 30% (Fig. 5, E and F). The longevity of crystal contact in wave-dominated flow is a consequence of individual crystals drifting along in the wave field, exhibiting only a slight wobble around a largely stable alignment (Fig. 5E). In contrast, crystals in pure shear at 5 to 15% crystal fraction tumble along, rotating even after coming in contact with another crystal (Fig. 5A).


Our models demonstrate that the olivine aggregates found in picritic scoria erupted at Kīlauea Iki in 1959 contain several subtle but distinct clues that allow us to infer specific quantitative attributes of the flow field at the time of aggregate formation. Here, we focus on understanding the following three observations: First, the aggregates are small with most containing only two to four similarly sized crystals (8, 10); second, an unexpectedly high percentage (≈37%) of aggregates exhibit a large misalignment angle between crystals (see Fig. 1E) (810); and third, the contact zones between the olivine crystals are partially or completely overgrown (see Fig. 1, A and B) (8, 10).

The first of these observations, the small aggregate size, suggests that aggregates formed at low crystal fraction. Both analog laboratory experiments of settling clay particles with shapes inspired by the Kīlauea Iki olivines (10) and a virtual reanalysis of these experiments (12) reproduced the formation of small aggregates through synneusis. However, neither laboratory experiments nor simulations were able to reproduce large misorientation angles within the aggregates or long-lived contact between crystals (11, 12). The second and third observations hence require a distinct physical process that explains both the observed misorientations and the long-lived contact.

Here, we test the hypothesis that exposure of the olivine crystals to a pronounced, steady wave field in the volcanic conduit could be the elusive mechanism reflected in the data. We show that for crystals with an axisymmetric symmetry like the Kīlauea Iki olivine there are two preferential alignment angles ±ϕ* with respect to the horizontal direction in the conduit. As a consequence, a wave-dominated flow would contain an approximately equal share of crystals aligned onto +ϕ* and −ϕ* angles. When neglecting the long-range, hydrodynamic interactions between crystals, aggregates would hence form either between crystals that are aligned onto the same angle, leading to an approximately zero misorientation angle, or between crystals that are aligned onto opposing angles, creating a misorientation angle of 2ϕ*. Since both alignment angles are equally likely, both cases should arise in about equal proportion. In the observational data, the first case would be represented by both the {110}-{110} and {010}-{010} contacts, since we do not resolve the crystallographic faces with sufficient detail to distinguish between the two, and the second case would be represented by the {021}-{021} contacts. Both classes of contact constitute roughly 40% of the total aggregates in the data, in consistency with model expectations.

While the analytical model is appealing for providing a process-based explanation for the formation of large misorientation angles, its neglect of the hydrodynamic interactions between crystals is difficult to justify when studying crystal aggregates. To quantify the conditions under which the analytically derived insights apply, we complement our analytical analysis with direct numerical simulations. Our simulations demonstrate that crystal-crystal interactions begin suppressing preferential alignment at average crystal fractions as low as a few percent (Fig. 6A) and all but eliminate preferential alignment at or above 10% (Fig. 6, B and C). As a consequence, we can only reproduce a substantial proportion of aggregates with large misorientation angles (10) for flows with small crystal fractions of approximately 5%. We emphasize that our modeling does not capture the geochemical processes creating overgrowth of olivines. Our simulations only focus on the fluid dynamical interactions between crystals and flow and hence provide snapshots of transient crystal contact, where individual crystals are able to separate from each other or drift together freely. As a consequence, our simulations tend to underestimate the total percentage of crystals clustered as compared to a scenario where crystals can overgrow.

Our simulations can, however, identify the necessary fluid-dynamical conditions for long-lived contact between crystals, at least at low crystal fraction. We find that the large proportion of aggregates with large misorientation angles is not the only observation that indicates a wave-dominated flow field in the volcanic conduit. Wave-dominated flow conditions in the conduit are also necessary to explain long-lived crystal contact, which is a prerequisite for the observed overgrown contact zones (Fig. 6, D and E). Contrary to wave-dominated flow, crystal contact in pure shear is short-lived because crystals continue to rotate, albeit at a reduced rate, after coming in contact. We emphasize that our finding of short-lived crystal contact in pure shear applies only at the relatively low crystal fractions considered here, when there is still enough space for crystals to rotate. As the crystal fraction increases, rotations are increasingly suppressed and long-lived crystal contact emerges in shear-dominated flow as well (12, 38). However, the crystal aggregates formed under these conditions tend to be larger in size, resembling force chains or fragments from networks (12, 38). Crystal aggregates formed in shear at intermediate crystal fraction are hence observationally distinct from the Kīlauea Iki aggregates, most of which contain two to four individual crystals.

Together, our modeling suggests that the Kīlauea Iki aggregates formed at low crystal fraction (5%) during exposure to a steady wave traveling downward into the conduit with wave-related shear exceeding the mean background shear by approximately a factor of three or more. In a volcanic conduit, the flow conditions that appear to have left their imprint on the Kīlauea Iki olivines arise naturally at the interface between gas-rich, buoyant, and degassed, heavy magma in bidirectional conduit flow (see Fig. 2). Crystals in the vicinity of this interface are exposed to a steady wave field with relatively slow background flow (22). The downward direction of the inferred wave is indicative of the portion of the conduit occupied by gas-rich, buoyant magma and its ascent speed (22, 24). The finding that misorientation angles are sensitive to the flow configuration highlights that tracking misorientation angles over time could help reveal if and how changes in conduit flow relate to eruptive behavior, even beyond the broad distinction between unidirectional or bidirectional flow we focus on here.

The presence of an interfacial wave not only produces large misorientations in crystal aggregates but also clears the interface of crystals (see Fig. 5, A to C), creating a locally low crystal fraction that is much smaller than in the rest of the conduit. Our simulations hence suggest that the crystal content within the conduit is spatially variable, even at a given depth section and under steady-flow conditions. More specifically, our analysis demonstrates that the samples analyzed here likely originated at the interface between the ascending and descending magmas and hence represent the low-end of observed crystal fractions. The lowest crystal fractions documented by Richter and Murata was indeed as low as 4%, but crystal contents overall span a wide range from 4% to 30%. However, this broad range is at least partially a consequence of the original samples representing all 17 episodes of the summit eruption during which both magma composition and crystal content changed notably (26, 39). The subset analyzed by Schwindinger (8) and used here are all scoria and mostly erupted during the relatively hot and crystal-poor first eruptive episode. It is hence conceivable that the Schwindinger samples would fall on the low range of the crystal fractions documented by Richter and Murata (26, 39).

Our modeling implies the testable prediction that samples from other portions of the Kīlauea Iki eruption, such as those erupted on the volcano flank or those probed in the summit lava lake (26, 39), do not contain an equally high percentage of aggregates with large misorientation angles as the samples (10) discussed here. We also posit that aggregates found in the lava flows created by the first fire-fountaining episode, which is best represented in the samples discussed here (10), would contain more crystals in total, but fewer crystals separated by large misorientation angles. A recent reanalysis of Kīlauea aggregates from multiple eruptions including but not limited to the 1959 data (9) shows that samples from other locations such as the East and Southwest Rift Zones contain much smaller percentages of aggregates with large misorientation angles (9), which is consistent with our model results. While encouraging, a more complete test of our interpretation would require comparing the petrography of crystal aggregates in scoria as opposed to lava flow samples from the same eruptive episode.

If our interpretation here is correct, then the Kīlauea Iki aggregates might be the first direct evidence of volcanic bidirectional flow. Originally proposed in the context of Kīlauea volcano (16), where backflow during eruptive activity including the 1959 eruption can be observed directly (40), bidirectional flow has been proposed to explain the long-term stability and endogenous growth of persistently active volcanoes. It is motivated by the vast quantities of gas escaping many persistently active volcanoes passively, which are difficult to reconcile with the exclusive degassing of magma extruded during eruptive and effusive episodes (17, 18, 4147). Since much more magma is degassed than erupted, magma must be recirculated creating bidirectional flow (1618, 20, 23, 48). However, while bidirectional flow is appealing for its ability to explain the large heat and gas fluxes characteristic of persistently active volcanoes (1618, 20, 23, 48), it does not offer an obvious explanation for why eruptions occur instead of all of the gas escaping passively.

One possibility for generating eruptions from bidirectional flow is through the ascent of a separate gas phase (23). The possible existence of an almost crystal-free interface around the uprising magma in the volcanic conduit, as inferred from our modeling, suggests another potential mechanism: Since strain localizes in crystal-poor regions, distributed, bidirectional flow throughout the conduit can transition to localized sliding at the crystal-poor interface (38). The transition entails a rapid increase in slip speed (49, 50), concentrating most of deformation on the sliding interface. The concentrated deformation favors gas nucleation (51) and further localization and degassing (52, 53), potentially to the degree that fragmentation occurs on the sliding interface (52, 53). The localized deformation and vesiculation on the sliding interface could also explain why magma from the interfacial environment is erupted as scoria. Over time, the slip-enabled, rapid ascent of magma creates an increasing mass imbalance in the volcanic conduit, implying that the fire fountain can only be sustained for a limited time.

This interpretation suggests that the 17 fire-fountaining episodes that modulated the 1959 summit eruption could have been created by processes that are not fundamentally different from the slip-enabled plug eruption assumed in prior conceptualizations of large-scale, eruptive activity (49, 50). However, there are several key differences. Instead of a solid plug, the magma ascending in the center of the conduit is merely more crystal-rich than the bidirectional flow interface. The small but finite difference in crystallinity is associated with a notable difference in effective viscosity (38), favoring slip at the bidirectional flow interface rather than at the conduit walls (49, 50). Also, the onset of slip is not driven by increasing pressure at depth (49, 50) but emerges from a dynamic slip instability on the bidirectional flow interface at relatively shallow depth. Needless to say, we can merely hypothesize the existence of this mechanism based on our analysis of the Kīlauea Iki aggregates. A more targeted model would be required to formally test the hypothesis posed here.

We also emphasize that our analysis here neither refutes nor confirms the idea that fire fountains might be created by the ascent of large gas pockets in the conduit (5456). Both types of basic conduit-flow models, unidirectional and bidirectional (see Fig. 2), are compatible with the ascent of a separate gas phase superimposed onto the background flow field. Nonetheless, the background flow field matters, because bidirectional flow entails an additional interface already prone to instability (22, 24, 25) and thus implies an inherently more dynamic view of volcanic conduits than unidirectional flow (23). The interactions with a separate gas phase could further amplify the existing interfacial instability.

While a separate gas phase almost certainly exists in the volcanic conduit (57), we do not consider it likely that the wave we infer here was generated by the passage of a large gas pocket. The main reason is that the inferred wave field must be steady over enough time such that the crystals may preferentially align and overgrow. It is not clear that a rapid disruption, such as the passage of a large gas volume (55) or a rockfall (7), could generate a steady enough wave field. Furthermore, the inferred wave field must have preceded the eruption, since most of the samples discussed here were collected during the first episode of fire fountaining and must have aggregated well before the 1959 eruptive sequence began.

Needless to say, more work is needed to fully test the inferences proposed here and the implications they might have for the eruptive mechanism generating fire-fountaining. Existing steady-state models merely demonstrate that a steady, downward-traveling wave is realistic (22, 24). It is not clear that conduit flow was at steady state before the onset of the 1959 eruption, and the assumption of constant fluid properties commonly made in the context of laboratory studies of core-annular flow (18, 21, 22) is inadequate for volcanic systems. We also note that we impose a relatively simple wave in our simulations to enable comparison with the linear waves assumed in our analytical model. The actual wave field in the conduit is likely more nonlinear, which may have more pronounced or additional effects on the crystals and the flow. In addition to further modeling work, systematically tracking crystal misorientation over time could enable us to identify changes in the conduit flow conditions and explore the degree to which these are linked to the waxing and waning of lava fountaining episodes.

Here, we leverage the subtle imprint left by volcanic conduit flow on olivine aggregates erupted in picritic socria at Kīlauea Iki in 1959 to infer several quantitative attributes of the pre-eruptive conditions in the conduit. We show that crystal aggregates with the observed characteristics form naturally in bidirectional, but not unidirectional flow. Even within the context of bidirectional flow, we show that specific flow conditions are required to reproduce the observations. We hence argue that the collected scoria samples are not representative of the conduit as a whole but are indicative only of the environment near the interface between gas-rich, buoyant, and degassed, heavy magma. The Kīlauea Iki samples might hence hold valuable clues for understanding the processes controlling fire-fountaining episodes, because they highlight that the bidirectional flow interface has low crystal fraction and is hence prone to transition from slow, distributed flow to rapid, localized sliding. In conclusion, our findings demonstrate that conduit flow information can be encoded on erupted crystal aggregates, and therefore, we suggest that tracking crystal misorientation throughout the different stages of future Kīlauea eruptions could shed light on conduit flow conditions over time.


Analytical model

The behavior of nonspherical particles, or crystals, has been studied under a number of idealized flow scenarios (30, 31). A small inertia-less particle in a fluid environment rotates according to Jeffery’s equation (32), which states that the particle’s angular velocity is a function of the flow’s rate of strain and rate of rotation and particle shape. The rate of rotation and strain can be thought of as competing effects, where the former tends to rotate the particles, and the latter tends to align the particles.

Under a steady, laminar shear flow, an ellipsoidal crystal in creeping flow will continually rotate undergoing what is referred to as a Jeffery orbit (32). The period of the orbit is a function of the ellipsoid’s aspect ratio λ and the flow’s shear rate γ. The angular velocity of the crystal is not constant, resulting in crystals aligning with the mean shear flow direction most of the time. A progressive wave superimposed on the shear flow is one way to break the limiting orbit.

The alignment of ellipsoidal particles has also been analyzed under traveling waves (13, 35). In waves described by linear wave theory, a small particle tends to a preferential alignment described by a stable limit cycle about a fixed point (13, 35). In contrast to the shear flow case, this preferential alignment is not necessarily parallel to the wave direction.

Examples of crystal alignment over time in pure-shear as compared to wave-driven flow are shown in Fig. 3 (A and B) in the main manuscript. The bidirectional conduit flow model is characterized by a sheared interface with a steady wave. It hence represents the case where both flow features, shear and wave, are superimposed. In contrast, the unidirectional flow model represents a steady shear flow. We derive an analytical model to describe the preferential alignment of nonspherical particles in a flow that combines a wave and a shear component and determine the theoretical limits of when crystals align preferentially.

Flow description. The flow field of interest is defined by the presence of a traveling wave train and a steady sheared current. We idealize the flow as a two-dimensional wave along an interface with a superimposed steady laminar shear flow, or Couette flow, as depicted in Fig. 2. We examine waves described by linear wave theory. This flow model is helpful because it has an analytical description for the velocity field generated by the waves and provides a theoretical framework to examine the alignment of nonspherical particles in response to fluid forcing.

Linear wave theory assumes small amplitude waves traveling in a two-dimensional, incompressible, and irrotational flow. In reality, the interface waves in the conduit are nonlinear (24, 58). We justify the use of linear wave theory in this analysis by noting that the wave-induced preferential alignment has been shown in nonlinear waves as well (59) and is confirmed by our direct numerical simulations that complement our analytical model and include viscous forces. The wave and shear system can be thought of as a flow system with both time-averaged mean and fluctuating strain and rotation components. We argue that the type of wave is not as important as the relationship between these mean and fluctuating flow components and that linear wave theory is a reasonable starting point.

We idealize the wave surface η as the deflection from the mean interface, x = 0 asη=Acos(kzωt)(2)where A is the wave amplitude, k is the wave number (associated with the wavelength L as k = 2π/L), and ω is the wave frequency (related to the wave period Tw as ω = 2π/Tw). The wave travels vertically in the positive z direction within our coordinate system. The horizontal coordinate is x such that the wave is aligned into the x-z plane. A negative wavenumber indicates a wave traveling downward in the negative z direction. Assuming small amplitude waves, i.e., kA ≪ 1, linear wave theory gives velocities induced by the wave around the interface as followsuw=Aωcosh[k(x+δ)]sinhkδcos (kzωt)(3)ww=Aωsinh[k(x+δ)]sinhkδsin(kzωt)(4)where uw = dz/dt is the vertical wave velocity, ww = dx/dt is the horizontal wave velocity, and δ is the width of the fluid between the wave interface and the wall. We will consider particles only on one side of the wave interface in this analysis; however, because of the symmetries of the flow, our results can be generalized to the entire conduit.

The shear flow idealized as a steady laminar shear flow has the following velocity fielduγ=U0+xγ(5)wγ=0(6)where uγ is the vertical velocity because of the shear flow, wγ is the horizontal velocity because of the shear flow, U0 is the velocity at the interface, and γ is the constant shear rate duγ/dx. In this case, we assume that the steady laminar shear flow can be superimposed on the velocity field induced by the waves in the frame of reference of U0. Thus, the total velocity field isu=Aωcosh[k(x+δ)]sinhkδcos(kzωt)+xγ(7)w=Aωsinh[k(x+δ)]sinhkδsin(kzωt)(8)

For simplicity, we assume that the shear component only affects the particle’s alignment but leaves its trajectory unchanged. Therefore, the particle’s relative trajectory is solely a function of the wave-induced flow and the mean shear flow. In other words, we assume that the difference in velocity in the mean flow across the particle’s orbit is much less than the velocities induced by the wave field. The wave orbit scales with wave amplitude A, so the velocity difference in the shear flow across the wave orbit will be γA; the wave velocity amplitude scales with ωA. Thus, the constraint on the flow is γA ≪ ωA, or more simplyγω(9)

Crystal description. We analyze the behavior of disperse crystals and neglect crystal-crystal interactions to isolate the flow-crystal interactions. We assume that the crystals are passive flow tracers, infinitesimal and inertia-less, so their particle Reynolds number Rep = 0, where Rep = ur, ur is the relative speed between the crystal and the fluid, d is the length scale of the crystal, and ν is the kinematic viscosity of the fluid. We also assume that the Stokes number St = 0, where St = τpf, τp is the relaxation time of the crystal, and τf is the shortest time scale of the flow. The relaxation time of the crystal represents the time it takes the crystal to react to the flow, so τp < < τf implies a crystal instantaneously reacts to the flow. Kīlauea Iki olivines are thought to be characterized by Rep ≪ 1 and St ≪ 1 (11).

Assuming the particles are spheroidal in shape as depicted in Fig. 3C, they will rotate with the flow according to Jeffery’s equation (32). Defining p as a unit vector that points along the particle’s axis of symmetry, Jeffery’s equation givesp·i=Ωijpj+ε(SijpjpipjSjkpk)(10)where Ωij is the flow’s rate of rotation tensor and Sij is the flow’s rate of strain tensor.

For simplicity, we have assumed the particles are axisymmetric. However, olivine crystals are typically non-axisymmetric particles that can be idealized as tri-axial ellipsoids. Tri-axial ellipsoids are defined by three unique radii lengths a, b, and c, where abc. The Kīlauea Iki olivines span a spectrum of average radii from prolate to oblate (8, 10) with a typical aspect ratio between 3.5:1 and 2:1 between the largest and smallest axis. Since previous work found that tri-axial particles still preferentially align under waves according to their most extreme aspect ratio (59), we limit our analysis to the simpler case of axis-symmetric particles.

The preferential alignment induced by waves is caused by the manner in which a particle samples the flow under a wave. The process is the angular analog to Stokes drift (13), where Stokes drift is the induced second-order Lagrangian transport of particles under waves with zero Eulerian mean flow (60). This preferential alignment is described by a steady limit cycle of the particle’s alignment angle ϕ, where ϕ is the alignment angle of the unit vector p, and ϕ=cos1px. The limit cycle is defined by the wave-averaged alignment, ϕ¯ where we have defined the wave-period average by, for exampleζ¯=12π02πζd(kzωt)(11)

The preferential alignment is defined as the fixed point of ϕ¯ where * denotes a fixed point. Under waves, it is solely a function of the particle’s eccentricity (59)ϕ¯*=12cos1ε(12)

Following a similar derivation to the pure wave case (13), we now include the effect of the shear flow on the particle’s alignment. Linear wave theory produces an irrotational flow, but the addition of a uniform shear flow makes Ωij ≠ 0. The nonzero rate of rotation terms now are Ωxz = γ/2 and Ωzx = − Ωxz.

The particle’s alignment is determined by Jeffery’s equation (32). In this coplanar wave and shear flow system, the Jeffery equations written in component form simplify top·x=εSxxpx(1(px2pz2))+εSxzpz(12px2)+Ωxzpz(13)p·y=pyε[Sxx(px2pz2)+2Sxzpxpz](14)p·z=εSxzpx(12pz2)εSxxpz(1+(px2pz2))Ωxzpx(15)

In the pure wave case, there is an unstable fixed point at py = 1, px = pz = 0. In this case we can derive the stability criterion and find that the eigenvalues Λ of the Jacobian about this fixed point are equal toΛ=±ε2(Sxx2+Sxz2)Ωxz2(16)

When Ωxz2<ε2(Sxx2+Sxz2), there always exists one real, negative eigenvalue, which implies that the fixed point is still unstable. On the other hand, when Ωxz2 is large and the argument of the square-root is negative, the eigenvalues of the Jacobian become purely imaginary. In that case, the fixed point is neutrally stable.

By numerically integrating these equations, we find that the particles do still settle onto a stable limit cycle under both waves and steady shear, unless the shear rate exceeds a critical threshold and Jeffery orbits reemerge. To determine the critical shear rate for this bifurcation to occur, we follow a similar derivation to that with only waves (13). To simplify the analysis, we convert our equations to a spherical coordinate system. When the particles are aligned within the x-z plane such that py = 0, the azimuthal angle θ = 0, and the alignment of the particle depend solely on the alignment angle ϕ=cos1px.

We obtain an equation that relates ϕ to the strain tensor terms, denoted by Sγ for shear flow components and Sw for wave components, and the flow’s rate of rotation Ωxzϕ·=ε(Sxxw+Sxxγ)sin 2ϕ+ε(Sxzw+Sxzγ)cos 2ϕ+Ωxz(17)

Given Sxzγ=Ωxz and Sxxγ=0, and using Eq. 2, we can plug in the expressions for the strain and rotation tensor termsϕ·=εkAωsinh kδ{sinh[k(x+δ)]cos(kzωt) cos 2ϕcosh[k(x+δ)]sin(kzωt)sin 2ϕ}+γ2(1+εcos 2ϕ)(18)

The full dynamical system is described by this governing equation for ϕ and the simplified relative equations of motion for the particle because of the wavesx·=Aωsinh [k(x+δ)]sinh kδsin (kzωt)(19)z·=Aωcosh [k(x+δ)]sinh kδcos (kzωt)(20)

Preferential alignment. The derived system of equations (18, 19, and 20) is nonlinear and nonautonomous and cannot be directly solved. To find a solution to this system, we use a Taylor expansion. We Taylor expand about xo and zo, which are the coordinates at the center of the wave orbital and follow the methodology described in (59). We include the leading-order finite-wave-amplitude correction to the crystal’s motion and alignment by using a Taylor expansion for Eq. 18. In this case, we expand about the crystal’s horizontal and vertical orbital motion (ξ and ζ respectively). For a function f, this expansion givesf(x,z,t)f(xo,zo,t)+ξfz|(xo,zo)+ζfx|(xo,zo)(21)whereξ(t)=Asinh [k(xo+δ)]sinh sin (kzoωt)(22)ζ(t)=Acosh [k(xo+δ)]sinh kδcos (kzoωt)(23)

Accounting for this first Taylor expansion term, we obtain an approximate expression for the crystal’s alignment as a function of tϕ·kAωεsinh kδ[cosh[k(xo+δ)]sin(kzoωt)sin2ϕ+sinh[k(xo+δ)]cos(kzoωt)cos2ϕ]+(kA)2ωε2sinh2kδ[sin[2(kzoωt)]sin2ϕ+sinh[2k(xo+δ)]cos2ϕ]+γ2(1+εcos2ϕ)(24)

Eq. 24 has one independent variable, t, but now includes a second order term, (kA)2. Because of the complexity of this equation, we implement a change of variables to enable analytical progress following the methods in (13). Let χ = sin 2ϕ and ψ = cos 2ϕ, and assume solutions for χ and ψ take the form of a fast mode that oscillates at the wave frequency ω and a slow mode that changes over a longer time scale, where B is some constantχ(t)=χ¯Bsin (kzoωt)(25)ψ(t)=ψ¯+Bcos (kzoωt)(26)

This change of variables appropriately recovers the dynamics of the full equations, as verified through a comparison of numerical integrations using the expressions for χ and ψ and the original ϕ equation. After changing variables, we average the equations over a wave period. This leaves us with an expression for ψ¯· as a function of the flow and particle parameters, and χ¯ and ψ¯ψ¯·=2kAωεsinh kδB(cosh(k(xo+δ))+12sinh(k(xo+δ)))χ¯(kA)2ωεsinh2kδsinh(2k(xo+δ))χ¯ψ¯γ(χ¯+εχ¯ψ¯)(27)

The preferential alignment of the particle is associated with the fixed point of this equation. To solve for it, we set the equation equal to zero, divide by χ¯, and solve for the fixed point, ψ¯*0=2kAωεsinh kδB(cosh(k(xo+δ))+12sinh(k(xo+δ)))(kA)2ωεsinh2kδsinh(2k(xo+δ))ψ¯*γ(1+εψ¯*)(28)

Rearranging, we findψ¯*=2kAωεsinh kδB(cosh(k(xo+δ))+12sinh(k(xo+δ)))γ(kA)2ωεsinh2kδsinh(2k(xo+δ))+εγ(29)

From DiBenedetto and Ouellette (13), the value of the coefficient B isB=kAε2sinh Kδ(sinh(2k(xo+δ))cosh(k(xo+δ))+12sinh(k(xo+δ)))(30)

After plugging this expression in for B in we are left withψ¯*=(kA)2ωε2sinh2kδsinh(2k(xo+δ))γ(kA)2ωεsinh2kδsinh(2k(xo+δ))+εγ(31)

In the case of zero shear, when γ = 0, we recover the purely wave result, ψ¯*=ε from (13). To simplify this new expression, we can use our physical intuition about the problem. Particle alignment is controlled by the flow’s velocity gradients. It has been shown that the preferential alignment is the angular analog to Stokes drift (13). Thus, it follows that the alignment would depend directly on the gradient of the Stokes drift velocity. The Stokes drift velocity ust to second order isust=kA2ωcosh(2k(xo+δ))2sinh2kδ(32)

The associated Stokes drift shear, γw = dust/dx is defined asγw=dustdx=(kA)2ωsinh2kδsinh(2k(xo+δ))(33)

If we plug γw into Eq. 31, we see that it simplifies, yieldingψ¯*=γwε2+γγwε+εγ(34)

Last, let 𝒟 = γ/γw, a ratio of the strength of the mean shear to the wave effects in this equation. Our solution simplifies toψ¯*=ε+D/ε1+D(35)

At this point, we can apply the change of variables, leaving the preferential alignment of the particle as a wave-period averaged polar angle that is solely a function of ε and 𝒟, resulting in Eq. 1.

Numerical model

To quantify the effect of crystal-crystal interactions on alignment, we compare our analytical results to direct numerical simulations based on the methodology derived by Qin et al. (12). The main added value of the direct numerical simulations is in capturing the long-range hydrodynamic interactions between crystals and model crystal collisions and aggregation. The numerical simulations also capture nonlinear wave effects, which we neglect in the analytical model.

Governing equation. Our numerical model solves for conservation of mass and momentum in two dimensions. We assume that the viscosity of the crystal-free magma is Newtonian and high enough for inertial effects to be negligible. We also assume that the magmatic liquid is incompressible. Therefore, the governing equations in the liquid phase are the incompressibility condition and the Stokes equation·v=0 and(36)p+η2v+ρg=0(37)where v represent the velocity, p the pressure, g the gravitational acceleration, ρ the density, η the dynamic viscosity of the crystal-free magmatic melt.

At the low to intermediate crystal fraction considered here, the crystals behave approximately like rigid bodiesMpdVpdt=Fp+Mpg(38)d(Ipωp)dt=Tp(39)dXpdt=Vp(40)where Mp represents the mass of the crystal, Vp the velocity of its center of mass, Xp the position of its center of mass, ωp its angular velocity, and Fp and Tp are the hydrodynamic force and torque exerted onto the crystal by the surrounding fluid. The angular moment of inertia tensor of the crystal is denoted by Ip.

Boundary conditions. In our simulations, we consider three different flow fields: a linear shear flow, a linear shear flow with a superimposed wave, and a pure-wave field. In pure shear, we define the shear speed, U0, to be zero on the left side of the domain, which represents the conduit wall. We compute the shear speed on the other three sides based on Eq. 5 to ensure that the shear-rate is the same as in the other two flow fields considered. For the pure wave field, we compute the vertical and horizontal speeds at the domain boundaries from Eq. 2 at x = 0 and x = − δ, to boundaries. We set the length of the computational domain to be same as the wavelength, L, and apply the solution of Eq. 2 at z = 0 and z = L to the upper and lower boundaries, respectively. To impose the same speed scale in pure shear as in the wave field, we enforce the same shear rate, γ, in both cases whereγ=duγdx=Aωcosh kδsinh kδδ(41)

When a traveling wave train is superimposed on pure shear, we compute the speed on the domain boundaries from Eq. 7. In all three flow fields, we allow the crystals to pass through the bottom and top boundaries periodically.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank P. Wallace for the suggestion to study the Kīlauea Iki olivines, K. Allison for contributions early on in the project and for producing several panels D and E in Fig. 1, G. Mahood, K. Cashman, and C. Culha for constructive edits of the manuscript text, and N. Ouellette for helpful discussions concerning the analytical model. We also thank E. Breard and one anonymous reviewer for valuable suggestions that helped improve the manuscript. Funding: M.D. acknowledges support the Stanford Gerald J. Lieberman Fellowship and the Postdoctoral Scholarship from Woods Hole Oceanographic Institution. Author contributions: M.D. derived the analytical model. Z.Q. performed the direct numerical simulations. J.S. conceptualized the study and advised the two model components. M.D. and J.S. wrote most of the text. All authors reviewed the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper. Open-source codes for reproducing the model results reported in this article are available from the GitLab repository of the SIGMA research group at Stanford University: Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article