Crystal and volatile controls on the mixing and mingling of magmas

The mixing and mingling of magmas of different compositions are important geological processes. They produce various distinctive textures and geochemical signals in both plutonic and volcanic rocks and have implications for eruption triggering. Both processes are widely studied, with prior work focusing on field and textural observations, geochemical analysis of samples, theoretical and numerical modelling, and experiments. However, despite the vast amount of existing literature, there remain numerous unresolved questions. In particular, how does the presence of crystals and exsolved volatiles control the dynamics of mixing and mingling? Furthermore, to what extent can this dependence be parameterised through the effect of crystallinity and vesicularity on bulk magma properties such as viscosity and density? In this contribution, we review the state of the art for models of mixing and mingling processes and how they have been informed by field, analytical, experimental and numerical investigations. We then show how analytical observations of mixed and mingled lavas from four volcanoes (Chaos Crags, Lassen Peak, Mt. Unzen and Soufrière Hills) have been used to infer a conceptual model for mixing and mingling dynamics in magma storage regions. Finally, we review recent advances in incorporating multi-phase effects in numerical modelling of mixing and mingling, and highlight the challenges associated with bringing together empirical conceptual models and theoretically-based numerical simulations.

modeling of mixing and mingling, and highlight the challenges associated with bringing together empirical conceptual models and theoretically based numerical simulations.

Introduction: Magma Mixing and Mingling and Volcanic Plumbing Systems
It is now widely accepted that magmas of different compositions can mix and mingle together (Blake et al., 1965;Eichelberger, 1980;Morgavi et al., 2019;Perugini & Poli, 2012;Snyder, 1997;Sparks & Marshall, 1986, Wiebe, 1987Wilcox, 1999). Textural consequences of mingling have long been observed (Judd, 1893;Phillips, 1880), although the earliest observations were not necessarily interpreted correctly (Wilcox, 1999), with heterogeneities interpreted as originating from metasomatism (Fenner, 1926) or solid-state diffusion (Nockolds, 1933). Advancements in geochemical analysis combined with an understanding of phase equilibria led to acknowledgment of mixing and mingling as key processes, alongside crystal fractionation, in producing the compositional diversity of igneous rocks (Vogel et al., 2008). In addition, interaction between magmas became recognized as a potential trigger for volcanic eruptions (Sparks et al., 1977). Evidently, understanding mixing and mingling processes is crucial for deciphering the evolution of igneous rocks and the eruptive dynamics of volcanoes.
Previous work has sometimes been flexible with regard to precise definitions of the terms "mixing" and "mingling." We here define mixing to be chemical interaction between two magmas that produces a composition intermediate between the original end-members (Bunsen, 1851). Chemical mixing proceeds by chemical diffusion (Lesher, 1994;Watson, 1982) and, if allowed to complete, leads to hybridization and homogeneous products . By contrast, mingling is the physical interaction of the two magmas, such as through convective stirring (e.g., Oldenburg et al., 1989) or chaotic where the bulk viscosities of the two magmas become closer, thereby facilitating mingling and mixing before continued crystallization of the mafic magma increases its viscosity.
Another scenario is mixing and mingling between partially molten silicic rocks and a hot, rhyolitic injection (Bindeman & Simakin, 2014), which is important for the formation of large, eruptible magma bodies containing crystals mixed from different portions of the same magma storage system (antecrysts; Bindeman & Melnik, 2016;Francalanci et al., 2011;Ubide et al., 2014a;Seitz et al., 2018;Stelten et al., 2015). In all cases, the physico-chemical changes and their associated timescales govern the style of mixing, the resultant textures, and the eruptive potential.
Evidence of mixing is preserved primarily at the microscale because the relatively slow rate of diffusion alone (Acosta-Vigil et al., 2012;Morgan et al., 2008) cannot redistribute chemical components over large spatial scales (Bindeman & Davis, 1999). Crystals, in particular, can preserve chemical records of changing storage conditions that can be associated with mixing. For instance, resorption zones and reverse zoning in plagioclase might indicate changes to more mafic melt compositions, possibly due to multiple mixing events (Hibbard, 1981;Lipman et al., 1997;Tsuchiyama, 1985). The mixing history can be determined by combining these observations with methodologies such as major-element (Rossi et al., 2019), trace-element , and isotopic analyses (Davidson et al., 2007), along with measurements from the bulk rock or other minerals. This can include timescales of mixing (Chamberlain et al., 2014;Rossi et al., 2019) and ascent , temperatures and pressures of mixing (Samaniego et al., 2011), and the relative contribution of processes such as fractional crystallization (Foley et al., 2012;Ruprecht et al., 2012;Scott et al., 2013).
Despite this, many studies continue to model mingling as taking place between two crystalfree fluids in a vat (Montagna et al., 2015). Such a picture is hard to reconcile with evidence from petrological analysis (Cooper, 2017;Druitt et al., 2012;Turner & Costa, 2007) and the lack of geophysical evidence for large extended bodies of melt (Farrell et al., 2014;Miller & Smith, 1999;Pritchard et al., 2018;Sinton & Detrick, 1992). It is therefore clear that the presence of crystals and volatiles, and their effect on magma rheology (Caricchi et al., 2007;Mader et al., 2013;Mueller et al., 2010;Pistone et al., 2012), must be accounted for when modeling mingling (Andrews & Manga, 2014;Laumonier et al., 2014).

Analogue Experiments
Early analogue experiments used non-magmatic fluids and particles to model magma mingling by injecting one viscous fluid into another (Campbell & Turner, 1986;Huppert et al., 1984Huppert et al., , 1986. These studies considered magmas as pure melts and demonstrated that large viscosity contrasts prohibit efficient mingling. Field observations that some mafic magmas became vesiculated in response to undercooling by the host magma (Bacon, 1986;Bacon & Metz, 1984;Eichelberger, 1980) motivated experiments focused on bubble transfer from one viscous layer into another, and demonstrated that the rise of bubble plumes could cause mingling (Phillips & Woods, 2001Thomas et al., 1993). Recent experiments have examined the effect of crystals on intrusion break-up. For example,  injected a particle-rich corn syrup (high density and viscosity) into a large, horizontallysheared body of particle-free corn syrup (low density and viscosity) to model the injection of cooling (partially crystallized) mafic magma into a convecting magma chamber. They found that low particle concentrations caused the injection to fragment and form "enclaves," whereas at high particle concentrations it remained intact and formed a coherent layer. These experiments further suggest that in the presence of a yield stress in the injected magma, the greater the bulk viscosity contrast the smaller the lengthscale of intrusion fragmentation, thus enhancing homogeneity at the macroscopic scale (Hodge & Jellinek, 2020). Although no analogue experiments have considered liquid injection into variably crystalline suspensions, experiments with gas injection into particle-liquid suspensions show a strong control of particle concentration and injection style, with a threshold between ductile and brittle behavior at random close packing (Oppenheimer et al., 2015;Spina et al., 2016).

High-Temperature and/or High-Pressure Experiments
Investigations of magma interactions in high-temperature and/or high-pressure experiments can be broadly divided into two categories. Static experiments consider the juxtaposition of heated magmas and study mixing resulting from the diffusion of different melt components Van der Laan & Wyllie, 1993;Watson & Jurewicz, 1984;Wylie et al., 1989). Fluid motion can still occur in these static experiments, as variable diffusion rates between elements can create density gradients that drive compositional convection (Bindeman & Davis, 1999). Additionally, because water diffuses much more rapidly than other components (Ni & Zhang, 2008), transfer of water from hydrous mafic magmas to silicic bodies lowers the liquidus temperature of the latter, leading to undercooling and the production of quenched margins in the mafic member, even without a temperature contrast (Pistone et al., 2016a). Bubbles that exsolve in a lower, mafic layer can also rise buoyantly into the upper layer, entraining a filament of mafic melt behind them (Wiesmaier et al., 2015). Such bubble-induced mingling can be highly efficient and has also been documented in natural samples (Wiesmaier et al., 2015). It has been proposed that a similar style of mingling can occur through crystal settling (Jarvis et al., 2019;Renggli et al., 2016).
Dynamic experiments apply shear across the interface between two magmas and reproduce mingling behavior. The shear can be applied in various ways, with a rotating parallel plate geometry (Kouchi & Sungawa, 1982, 1985Laumonier et al., 2014Laumonier et al., , 2015, a Taylor-Couette configuration (De Campos et al., 2004, 2008Perugini et al., 2008;Zimanowski et al., 2004), a Journal Bearing System (De Campos et al., 2011;, or by using a centrifuge . These experiments have produced a variety of textures from homogenous mixed zones to banding. When pure melts are used, the combination of diffusional fractionation and chaotic advection can produce phenomena such as doublediffusive convection  and reproduce nonlinear mixing trends for various major and trace elements (De Campos et al., 2011;Perugini et al., 2008).
Experimental results also suggest new quantities to describe the completeness of mixing, such as the concentration variance  and the Shannon entropy . Where crystals are considered, the presence of phenocrysts can enhance mingling by creating local velocity gradients and disturbing the melt interface (De Campos et al., 2004;Kouchi & Sunagawa, 1982, 1985). In contrast, other studies (Laumonier et al., 2014(Laumonier et al., , 2015 have shown that the presence of a crystal framework in the mafic member prevents mingling, whereas the presence of water can enhance mingling by lowering the liquidus temperature, and thus the crystallinity, of the magma (Laumonier et al., 2015). Sparks and Marshall (1986) developed the first simple model to describe viscosity changes caused by thermal equilibration of a hot mafic magma and a cooler silicic magma, and the resulting (limited) time window in which mingling/mixing can occur. More sophisticated models have simulated mingling between melts driven by double-diffusive convection (Oldenburg et al., 1989), compositional melting (Cardoso & Woods, 1996;Kerr, 1994), and the Rayleigh-Taylor instability (Semenov & Polyansky, 2017). Another group of studies has used single-phase models to simulate elemental diffusion and advection in a chaotic flow field (Perugini & Poli, 2004;Petrelli et al., 2006). These models reproduce naturally observed geochemical mixing relationships, including linear-mixing trends between elements with similar diffusion coefficients and large degrees of scatter when diffusion coefficients differ (Nakamura & Kushiro, 1998;Perugini & Poli, 2004). Interestingly, the simulations produce both regular and chaotic regions, which are unmixed and well mixed, respectively, and have been interpreted to correspond to enclaves and host rock (Petrelli et al., 2006). This framework has been extended to account for a solid crystal phase (Petrelli et al., 2016) by including a Hershel-Buckley shape-dependent rheology (Mader et al., 2013) and a parameterization of the relationship between temperature and crystallinity (Nandedekar et al., 2014). This body of work has demonstrated that chaotic advection can speed up homogenization.

Numerical Models
Models of mixing and mingling that consider two-phase magmas containing either solid crystals or exsolved volatiles often assume coupling between the phases. In this way, the solid or volatile phase can be represented as a continuous scalar field, and the resultant effect on rheology is accounted for through a constitutive relationship. For example, Thomas and Tait (1997) used such a framework to show that volatile exsolution in an underplating mafic magma could create a foam at the interface with an overlying silicic magma. Depending on the exsolved gas volume fraction and melt viscosity ratio, mixing and mingling could then proceed through foam destabilization, enclave formation, or a total overturn of the system. Folch and Martí (1998) showed analytically that such exsolution could lead to overpressures capable of causing volcanic eruptions. Recent finite-element models show that injection of a volatile-rich mafic magma into a silicic host can cause intimate mingling when viscosities and viscosity contrasts are low (Montagna et al., 2015;Morgavi et al., 2019). The combination of reduced density in the chamber and the compressibility of volatiles can (non-intuitively) lead to depressurization in the chamber (Papale et al., 2017), which is important for interpretation of ground deformation signals (McCormick Kilbride et al., 2016).
The effect of crystals on mixing and mingling has also been modeled by treating the crystals as a continuous scalar field. Examples include simulations of mixing across a vertical interface between a crystal suspension (30% volume fraction) and a lighter, crystal-free magma (Bergantz, 2000), and injection of a mafic magma into a silicic host with associated melting and crystallization (Schubert et al., 2013). The role of crystal frameworks in both the intruding and host magma is addressed by Andrews and Manga (2014), who model the role of thermal convection in the host, and associated shear stress on the intruding dike. If convection occurs while the dike is still ductile, then mingling will produce banding.
Otherwise, the dike will fracture to form enclaves. Woods and Stock (2019) have also coupled thermodynamic and fluid modeling to simulate injection, melting, and crystallization in a sill-like geometry.
Finally, isothermal computational fluid dynamic simulations have been used to examine the case of aphyric magma injecting into a basaltic mush. For sufficiently slow injection rates, the new melt percolates through the porous mush framework, whereas for faster injections, fault-like surfaces delimit a "mixing bowl" within which the crystals fluidize and energetic mixing takes place (Bergantz et al., 2015Carrara et al., 2020;McIntire et al., 2019;Schleicher et al., 2016). By explicitly modeling the particles with a Lagrangian scheme, it is possible to account for particle-scale effects, including lubrication forces (Carrara et al., 2019), that are neglected when using constitutive relations from suspension rheology. These simulations suggest that mushes with ≤60% crystals can be mobilized by injection, but neglect welded crystals or recrystallization of crystal contacts. Furthermore, geophysical observations suggest that mushes spend the majority of their lifetimes with much higher crystallinities (80%-90%; Farrell et al., 2014;Pritchard et al., 2018;Sinton & Detrick, 1992). Despite these limitations, recent simulations using the model have shown that the contrast between the intruding and resident melt densities, rather than bulk densities controls the morphology of intrusion (Carrara et al., 2020).

Chaos Crags
Chaos Crags comprises a series of enclave-bearing rhyodacite lava domes that erupted between 1125 and 1060 years ago (Clynne, 1990). The host lavas are crystal-rich, containing phenocrysts of plagioclase, hornblende, biotite, and quartz, whereas the enclaves are basaltic andesite to andesite with occasional olivine, clinopyroxene, and plagioclase phenocrysts in a groundmass of amphibole and plagioclase microphenocrysts (Heiken & Eichelberger, 1980).
Many, but not all, enclaves have fine-grained and crenulated margins, and all contain resorbed phenocrysts captured from the host (Figure 4a). Some phenocrysts in the host also show resorption textures (Tepley et al., 1999).

Enclave Groundmass Textures
The enclaves from all four volcanoes show both similar and contrasting textural features. At Chaos Crags, most enclaves have fine-grained and crenulate margins (Figure 4a; Tepley et al., 1999), although those erupted in later domes are more angular and lack fine-grained margins. Enclaves in Lassen Peak samples are subrounded to subangular with an equigranular texture (Figure 4b; Clynne, 1999). Many enclaves from the 1991-1995 eruption at Mt. Unzen have crenulate and fine-grained margins (Browne et al., 2006a), although some have angular edges and a uniform crystal size (Figure 4c; Fomin & Plechov, 2012). Similar features are observed at Soufrière Hills, with many inclusions being ellipsoidal ( Figure 4d) and some angular; most, but not all, have fine-grained, crenulate margins (Murphy et al., 2000). Both the size and volume fraction of enclaves increased during the eruption Plail et al., 2014Plail et al., , 2018.
In all localities, fine-grained margins and crenulate contacts are attributed to undercooling of the mafic magma due to juxtaposition against the much cooler felsic host (Eichelberger, 1980) and associated rapid crystallization of the mafic melt near the contact with the felsic host. These crystalline rims have a greater rigidity than the lower-crystallinity enclave interiors so that as the enclave continues to cool and contract, the rims deform to a crenulate shape that preserves the original surface area (Blundy & Sparks, 1992). Enclaves not exhibiting such quench textures are also found at all localities.

Plagioclase
The composition and texture of plagioclase crystals are extremely good recorders of magmatic processes because (a) their stability field in pressure-temperature-composition (P-T-X) space is very large in volcanic systems, and (b) compositional zoning modulated by changes in the P-T-X space is well preserved due to the relatively slow diffusion in the coupled substitution between Na-Si and Ca-Al (Berlo et al., 2007;Grove et al., 1984;Morse, 1984).
Texturally, plagioclase phenocrysts in the host lavas at all four localities comprise a population of unreacted, oscillatory zoned crystals with a smaller amount of reacted crystals that have sieved cores and/or resorption rims (Figure 5a; Browne et al., 2006b;Clynne 1999;Murphy et al., 2000;Tepley et al., 1999). Associated enclaves contain plagioclase xenocrysts incorporated from the host with sieved-texture resorption zones that consist of patchy anorthite-rich plagioclase and inclusions of glass (quenched melt). These reacted zones can penetrate to the cores of smaller crystals (Figures 5b,c), but in larger xenocrysts appear as a resorption mantle surrounding a preserved oscillatory zoned core ( Figure 5d). All xenocrysts are surrounded by a clean rim that is of the same composition as the plagioclase microphenocrysts in the enclave groundmass.

Interpretation of Textures and Chemistries
The common textural and chemical features of these volcanic systems suggest commonalities in the mixing and mingling processes. First, because enclaves from all volcanoes contain xenocrysts that originated in the host magmas, the mafic component must have been sufficiently ductile to incorporate these crystals during mixing. Plagioclase xenocrysts contain rounded, patchy zones with a sieved texture showing that both partial and simple dissolution occurred (Cashman & Blundy, 2013;Nakamura & Shimakita, 1998;Tsuchiyama, 1985), suggesting that the enclave magmas were undersaturated in plagioclase at the time of incorporation. Because up to 70% of the enclave groundmass consists of plagioclase microphenocrysts, this implies the mafic magmas were crystal-poor at the time of xenocryst incorporation.
Compositional variations of FeO and An in the plagioclase crystals provide further information on the relative compositions of the host and enclave melt at Soufrière Hills    (Ruprecht & Wörner, 2007). At Mt. Unzen, enclave microphenocryst and xenocryst rims show a strong positive correlation for the whole An range, whereas these phases at Soufrière Hills show a negative correlation for An > 75 mol% ( Figure 6). This difference is attributed to the absence of Fe-Ti oxide as an early crystallizing phase in the Soufrière Hills mafic end-member, which would cause FeO to increase in the residual melt as other phases precipitated until the point of oxide saturation . The lack of this inflection in the Mt. Unzen sample suggests that Fe-Ti oxides were present in the mafic magma prior to mixing, as suggested for the 1991-1995 eruption (Botcharnikov et al., 2008;Holtz et al., 2005).
Whereas the observed enrichment in FeO in enclave microphenocrysts, sieved zones in phenocrysts and xenocrysts, and xenocryst rims is likely due to crystallization from a more mafic melt, it is also possible that growth of these regions may be sufficiently fast for kinetic effects to play a role; if growth is faster than diffusion of FeO in the melt, then an FeO-rich boundary layer may develop around the crystals (Bacon, 1989;Bottinga et al., 1966;Mollo et al., 2011) that could also explain the enrichment. However, such a process would generate a negative correlation between FeO and An , not the positive correlation observed at Unzen and Soufrière Hills.
The contrasting textures of quartz in the host and enclaves also provide insight into the mingling/mixing process. Rounding of quartz xenocrysts, together with glass-filled embayments, suggests dissolution of quartz in the host. Conversely, quartz reaction rims comprising hornblende microphenocrysts, glass, and vesicles in the enclaves (Figures 3d, 7b) suggest that the dissolution-induced increase in the silica content (and H2O solubility) of the surrounding melt caused diffusion of H2O toward the quartz (Pistone et al., 2016a).
Whereas the presence of resorbed xenocrysts in enclaves suggests that there was time for crystals to be incorporated, and to react, before the enclave started to crystallize, the presence of fine-grained rims on some enclaves Browne et al, 2006a;Murphy et al., 2000;Plail et al., 2014;Tepley et al., 1999) implies rapid cooling and crystallization (chilling) of the mafic magma against the cooler silicic host (Bacon, 1986). Xenocrysts must therefore have been incorporated prior to the formation of the chilled margin, providing a limited temporal window for crystal transfer. A comparison of the thickness of xenocryst resorption zones at Mt. Unzen (Browne et al. 2006a) with those produced experimentally (Nakamura & Shimakita, 1998;Tsuchiyama & Takahasi, 1983;Tshuchiyama, 1985) suggests resorption on a timescale of days; this contrasts with thermal modeling (Carslaw & Jaeger, 1959) suggesting that enclaves should thermally equilibrate on a timescale of hours. Again, this requires incorporation of xenocrysts prior to intrusion disaggregation and enclave formation (Browne et al., 2006a). As all the considered volcanic lavas contain similarly resorbed plagioclase xenocrysts within enclaves of comparable sizes, it seems likely that this temporal constraint on the sequence of crystal transfer prior to enclave formation is generally true for the systems presented here.
Importantly, all locations also contain enclaves with unquenched margins Tepley et al., 1999) and equigranular textures (Browne et al., 2006a;Heiken & Eichelberger, 1980). Equigranular enclaves at Mt. Unzen have been interpreted as originating from disaggregation of the interior of the intruding magma, which cooled more slowly than the intrusion margin where porphyritic enclaves (xenocrysts-bearing, chilled margin) formed.
Similarly, at Soufrière Hills, the quenched enclaves may form from an injected plume of mafic magma, whereas unquenched and more hybridized enclaves form from disturbance of a hybrid layer at the felsic-mafic interface . Angular enclaves with unquenched margins may record the break-up of larger enclaves (Clynne, 1999;Fomin & Plechov, 2012;Murphy et al, 2000;Plail et al., 2014), which can return resorbed host-derived crystals to the host; this explains the presence of resorption zones in crystals in the host lavas Further support for enclave fragmentation comes from microlites that are chemically indistinguishable from enclave phases at Soufrière Hills . A possible method to determine whether equigranular enclaves form from a hybrid layer or disaggregation of larger enclaves is to examine the mineralogy of the crystals in the enclave.
The two different mechanisms will produce different degrees of undercooling within the enclave magma, which, in the hybrid-layer model, will depend on the relative proportions of the end-member magmas, and thus can produce different crystal assemblages/textures .

Conceptual Model of Magma Mixing and Mingling
The common features of the eruptive products described above suggest common aspects of mixing and mingling. Xenocrystic mafic enclaves with chilled margins, in particular, require that magma injection be accompanied by crystal incorporation from the host magma, as also suggested by a comparison of thermal timescales with the times needed to generate the observed thicknesses of resorption zones (Browne et al., 2006a). These constraints on the sequence of mixing processes have led to a similar conceptual model of mixing and mingling ( Figure 8; Browne et al., 2006a;Clynne, 1999;Murphy et al., 2000;Plail et al., 2014;Tepley et al., 1999) in which the mafic magma is injected as a fountain (Clynne, 1999) or collapsing plume  before ponding at the base of the silicic host (Figure 8a). Shear caused by the injection disrupts the interface between the two magmas, leading to the formation of blobs of hybridized magma with incorporated host crystals that then rapidly chill against the silicic host, preventing further hybridization Tepley et al., 1999). Heating of the host, in turn, causes partial melting, reducing the crystallinity and causing convective motions that disperse the enclaves. Meanwhile, at the mafic-silicic contact, a hybrid interface layer forms ( Figure 8b). As this layer crystallizes, second boiling drives fluid saturation; exsolved buoyant fluids produce a low-density, gravitationally unstable, interface layer that breaks up to form further enclaves (Figure 8c; Browne et al., 2006a;Clynne, 1999). As cooling propagates downward through the mafic body, enclaves can come from deeper portions resulting in more equigranular enclaves that lack chilled margins or xenocrysts (Brown et al., 2006a;Plail et al., 2014). Enclaves, once formed, can disaggregate. Disaggregation is shown by the presence of broken enclaves (Clynne, 1999;Fomin & Plechov, 2012;Tepley et al., 1999), host phenocrysts with resorption zones and Fe enrichment caused by previous engulfment in mafic magma (Browne et al., 2006b;Clynne, 1999;Humphreys et al., 2009;Tepley et al., 1999), and small clusters of enclave-derived microlite material within the host lavas .
Disaggregation allows for subsequent mixing of a type precluded during initial enclave formation, but the timing of disaggregation is poorly constrained. It could occur during highshear conditions in the conduit ; alternatively, disaggregation may be part of a continuous cycle of injection, enclave formation, and fragmentation ( Figure 8d) that gives rise to a continuously convecting magma storage region, which is sometimes sampled during a volcanic eruption (Browne et al., 2006a). Regardless, the dispersion of mafic groundmass into the host has implications for interpreting end-member compositions from petrologic studies Martel et al., 2006). Importantly, neglecting such transfer can lead to an underestimate of the initial silica content of the felsic member.

Quantitative Modeling of Crystal and Volatile Controls on Mixing and Mingling
Many conceptual models of magma mixing (e.g., Figure 8) have been produced based on petrologic evidence. However, quantitative models of magma mixing are limited. As described in Section 2.4, Sparks and Marshall (1986) first developed a simple model describing how thermal equilibration of a juxtaposed mafic and silicic magma led to rapid viscosity changes that inhibited mixing after a short time. Since then, models developed to account for the role of either crystals or exsolved volatiles have produced significant insights into mingling and mixing dynamics, but have failed to incorporate petrological data within quantitative frameworks. Here, we examine three models: Andrews and Manga (2014), who use continuum modeling and suspension rheology to model mingling resulting from dike injection into a silicic host; Bergantz et al. (2015), who model the injection of melt into a basaltic mush, resolving both fluid and granular behavior; and Montagna et al. (2015), who simulate the effect of exsolved volatiles on mafic injection. We compare the model assumptions and results, as well as their implications for interpreting petrological data.

The Model of Andrews and Manga (2014)
The model considers the instantaneous injection of a mafic dike into a silicic host, with a prescribed initial composition and temperature, and numerically solves the 1D heat equation.
Changes in the crystallinity and bulk viscosity of magmas with time are calculated using MELTS simulations (Asimow & Ghiorso, 1998;Ghiorso & Sack, 1995;) and viscosity models for melt (Giordano et al., 2008) and crystal-bearing suspensions (Einstein, 1906;Roscoe, 1952). If the viscosity of the host immediately juxtaposed with the dike decreases sufficiently, then the host starts to convect (as determined by a Rayleigh number criterion), which exerts a shear stress on the dike. If this shear stress exceeds the yield stress of the dike (which depends on its crystal content), the dike deforms in a ductile fashion and the model predicts banded products. Alternatively, if the yield stress exceeds the shear stress, then the dike fractures in a brittle fashion and enclaves form.
In this model context, the principal control on mingling dynamics is the development of crystal frameworks within the dike. Dike crystallization, in turn, is controlled by composition and temperature contrasts. For example, injection of hot, large, and wet dikes causes the silicic host to convect before a crystal framework forms in the dike. The resultant shear causes ductile disruption of the dike and intimate mingling of the two magmas, producing banding and, with time, homogenization. Small and dry dikes, by contrast, experience extensive crystallization before the host starts to convect and thus fracture to form enclaves.
The precise initial conditions (temperature, dike size, and water content) that determine mingling style are sensitive to the parameterizations used (e.g., critical Rayleigh number for convection), but the qualitative results are useful.
The principal limitation of the model of Andrews and Manga (2014) is that it assumes an instantaneous injection of the mafic dike and therefore neglects any mixing/mingling that occurs during injection itself. Instead, the dike is disrupted only by shear due to convection in the host. Indeed, the relative importance of shear due to injection versus shear due to convection remains a considerable unknown. The assumption that brittle fragmentation of the dike produces enclaves is supported by three-dimensional tomographic observations of enclaves from Chaos Crags, which have crystal frameworks that are lacking in banded pumices from Lassen Peak (Andrews & Manga, 2014). The inference is that these crystal frameworks created a yield stress such that the enclaves formed by solid-like fracturing and banded pumice by ductile deformation. However, this is in direct contradiction with the conceptual model presented above (Figure 8), which is based on field and petrographic observations that suggest enclaves form from fluid-like deformation of the mafic magma.
This contradiction highlights the extent to which conditions of enclave formation are unknown. Bergantz et al. (2015) The discrete-element model, which resolves both fluid and granular physics, considers the injection of a crystal-free magma into the base of a crystal mush at random loose packing (approximately 60% crystallinity). The response of the mush is governed by stress chains formed by crystal-crystal contacts. For sufficiently slow injections, the new melt permeates through the mush, which behaves as a porous medium. Once the injection speed is large enough to disrupt the stress chains, however, part of the mush can become fluidized to form a mixing cavity, which is an isolated region where the host melt, crystals, and new melt undergo overturning. The new melt then escapes from the cavity through porous flow into the rest of the mush. For still faster flow speeds, the stress chains orientate to create two fault-like surfaces at angles of about 60° to the horizontal that bound a fluidized region of the mush, within which extensive circulation occurs. Recently, this model has been extended to investigate the effect of a density contrast between the intruding and resident melts on the style of mingling (Carrara et al., 2020), showing that the intrusion geometry is controlled to first order by the contrast between the melt densities rather than the bulk densities.

The Model of
Although this model captures granular and fluid dynamics on the crystal scale and demonstrates the impact of varying the injection velocity, there are numerous outstanding questions. First, varying the crystallinity of the mush has not been addressed and will presumably affect the values of the injection velocity at which transitions between mingling styles occur. Furthermore, temporal and spatial variations in temperature (due to heat transfer or latent heat release), and therefore in viscosity and crystallinity, have not been considered.
Cooling and crystallization of the new melt should control the dynamics of the system, as will associated latent heat release. Finally, the geometry of the modeled magma reservoir (laterally homogenous layers) will affect the specifics of the mixing process, such as the orientation of the bounding faults, and it is not yet clear if the model scales to natural systems.

The Model of Montagna et al. (2015)
The two-dimensional finite-element model considers two vertically separated magma chambers that are superliquidus and connected by a narrow conduit. The upper chamber initially contains a felsic phonolite, and the lower chamber and conduit are filled with a mafic shoshonite, compositions chosen to represent eruptions from Campi Flegrei. H2O and CO2 exsolve as functions of temperature and pressure (Papale et al., 2006), whereas the transport of exsolved volatiles is modeled as a continuum scalar field satisfying a transport equation.
Bubbles are assumed to be sufficiently small that they are undeformable, and an empirical law is used to parameterize their effect on bulk viscosity (Ishii & Zuber, 1979). The shoshonite initially contains exsolved volatiles and so is lighter than the phonolite, creating an unstable density interface at the inlet to the upper chamber.
Upon initiation, a Rayleigh-Taylor instability develops at the inlet to the upper chamber, and a plume of light material rises into the chamber while the conduit is filled with a mixed, hybrid magma. Intimate mingling within the chamber is reminiscent of that created by chaotic advection (Perugini & Poli, 2004). The magma entering the upper chamber is a partial hybrid, and the pure parent shoshonite never enters the upper conduit. Intensive mingling occurs on a timescale of hours, promoted by a large initial density contrast and horizontally elongated chambers. Importantly, the reduction in density of the upper chamber can cause depressurization, which has implications for interpreting ground deformation signals (Papale et al., 2017).
Although an obvious limitation of the model is the two-dimensional domain, it seems reasonable that the results can be extrapolated to three-dimensional systems. A greater limitation is the restricted range of compositions and temperatures for which the model is valid. The end-member compositions are similar and superliquidus, so that both the absolute bulk viscosities (<3500 Pa s) and their contrast (factor of 7) are relatively low. This allows rapid mingling and entirely ignores the effect of crystals on the flow dynamics.

Comparison and Common Limitations
Both Andrews and Manga (2014) and Bergantz et al. (2015) focused on the effect of crystals, but a key difference in the two models is the initial condition. Andrews and Manga (2014) assume the instantaneous injection of a dike into an initial rheologically locked host, whereas Bergantz et al. (2015) simulate the flow of new melt into a melt-crystal mixture; they show that new melt flows permeably through a rheologically locked mush. The conditions that spatially constrain a mafic injection (e.g., as a dike) have not been defined. The two models also simulate the role of crystals differently. Andrews and Manga (2014) calculate the crystallinity of a magma at a given temperature and assume the presence of a crystal framework (and yield stress) above a threshold value. Bergantz et al. (2015) allow the crystals to form force chains through which stresses are transmitted , but they consider the system to be isothermal such that no crystallization occurs, a key feature of Andrews and Manga (2014).
Both models are limited in addressing the role of volatiles. Diffusion of volatiles from the mafic to felsic member can strongly influence the crystal composition and textures of the silicic member (Pistone et al., 2016a), whereas exsolution of volatiles leads to a reduction in bulk density that can drive convective motions in the mixing dynamics (Eichelberger, 1980;Montagna et al., 2015;Phillips & Woods, 2001;Thomas et al., 1993;Wiesmaier et al., 2015). The presence of exsolved volatiles also affects the magma rheology and requires the use of three-phase rheological models (Mader at al., 2013;Pistone et al., 2016b). One strategy is to treat the exsolved phase as a continuum scalar field and use a suspension model for bulk rheology (Montagna et al., 2015). However, as has been shown for solid phases (Carrara et al., 2019), small-scale effects can be overlooked by this approach, and explicit modeling of such phases may be needed to accurately constrain mixing/mingling processes.
Additional complications arise in the number of parameters required for a given model. For example, the Andrews and Manga (2014) model requires values for a maximum crystal packing fraction and a critical Rayleigh number for convection in the host. Constraining these parameters will require extensive experimental efforts involving both high-temperature/highpressure and analogue experiments.

Conclusions and Outlook for Future Research
We have reviewed progress in understanding magma mixing and mingling, focusing on volatile and crystal controls on mingling processes. Although field and petrologic observations of mixed and mingled products are numerous, models of these processes do not yet include the full range of observed complexities. In particular, conceptual models derived from observations (Browne et al., 2006a;Clynne, 1999;Plail et al., 2014;Tepley et al., 1999;) suggest very different dynamics to those from numerical models (Andrews & Manga, 2014;Bergantz et al., 2015;Montagna et al., 2015). To resolve this discrepancy, several key questions need to be addressed: 1. How do mixing and mingling occur within the framework of crystal mushes, and how does the volume fraction of crystals control the interaction dynamics? 2. How do volatiles, both exsolved and dissolved, affect mixing and mingling? What is the relative importance of chemical quenching (due to volatile diffusion) versus thermal quenching (due to heat diffusion)?
3. How much mingling/mixing takes place during intrusion of the mafic magma compared to that driven by later processes such as convection in the host or the buoyant rise of vesicular mafic/hybrid magma? Only by combining field and analytical observations with experimental (analogue and natural materials) and numerical modeling can we start to address these challenges.