Plumes of light magma coming from depth keep entering the shallow reservoir as discrete filaments, following irregular trajectories and showing typical convective patterns. The lighter material tends to rise into the chamber, thereby decreasing more and more its density as volatiles exsolve in lower-pressure environments; on the other hand, the denser magmatic mixture initially residing in the chamber sinks into the feeder dyke, increasing its density by the reverse process of volatile dissolution at higher pressures. The plumes thus progressively increase their buoyancy, enhancing their expansion and acceleration. During the rise, vortexes form at the head of the plumes and subsequent plumes interact among themselves, further favouring mixing. The dynamics creates complicated patterns that maximize the interaction among the two different magmatic mixtures (Petrelli et al. 2011). Mingling is evident for all simulated conditions both within the chamber itself and even more in the feeding dyke (Fig. 2), and it is strongly intensified by the chaotic patterns that form as a consequence of deep magma injection.

Independently from system geometry or density contrast at the interface, mingling is very efficient in the feeding dyke, more than inside the upper chamber. Figure 2 shows that since the very beginning of the simulations, the magma entering the chamber is already a mixture of the

two initial end-members, and not the pure shoshonitic composition.

As the dynamics proceeds, faster for higher density contrasts and sill-like setups, the gas-rich mixture tends to accumulate at the top of the chamber, thereby originating a stable density stratification that has indeed been testified at various magmatic systems (Arienzo et al. 2009). The stratification is more prominent in vertically elongated, dyke-like reservoirs (Fig. 2). The density profile along the vertical direction, evaluated averaging along horizontal planes (Fig. 3), illustrates that a quasi-stable profile is reached after some hours of simulated time.

As time proceeds, convection slows down due to smaller buoyancy of the incoming already mixed component, and the instability proceeds in time asymptotically: the more the two end-members have mingled, the less intense is convection.

The evolution of pressure in the system is highly heterogeneous in space and time. Alternating phases dominated by buoyancy and sinking at chamber inlet result in pressure fluctuations with periods of hundreds of seconds and amplitudes decreasing with time (Fig. 4). Typically pressure variations are smaller than 1 MPa; under these conditions, it is unlikely that rejuvenation can trigger eruption, as the stresses needed to create a pathway to the surface in the host rock are typically larger than that (Gudmundsson 2006).

<figure><img src="image_41.png" /><figcaption>**Fig. 3** Total mixture density averaged over horizontal planes as function of depth for simulation 1, at different times. The inset shows the upper 5 km of the domain; the</figcaption></figure>

black line represents the quasi-equilibrium density profile at the end of the simulation
<figure><img src="image_42.png" /><figcaption>**Fig. 4** Pressure variations as a function of time at a point on the boundary of the upper chamber, for simulation #1. The *upper* diagram shows the difference between the local</figcaption></figure>

pressure at current time and at time zero, while the *bottom* diagram shows the same quantity after subtraction of a detrending function (*red curve* in the *upper diagram*)

## 3.3 Ground Deformation

Determining the time-space-dependent ground displacement requires modeling the magma-rocks boundary conditions and the mechanical response of rocks, the latter depending on heterogeneous rock properties, presence and distribution of faults, interfaces, fluids, and volcano topography (e.g., O'Brien and Bean 2004). A first-order analysis performed here assumes magma-rock one-way coupling and adopts the Green's functions formulation for a homogeneous, infinite medium (Aki and Richards 2002).

We consider as point sources the fluid dynamics computational grid nodes located at the reservoir walls. As source time functions, we use the respective temporal evolutions of magmatic forces computed from pressures and stresses provided at those nodes by the numerical simulations of magma convection and mixing dynamics. Ground displacement at a series of virtual receivers is finally obtained by integrat-

ing, over all sources, the Green's functions associated with individual sources.

Continuity of pressure and stress is taken as the boundary condition along the non moving magma-rock interface. Physical properties of rocks are homogeneous averages that describe the volcanic edifices within the range of considered depths (<10 km, v<sub>P</sub> = 3000 m/s; v<sub>P</sub>/v<sub>S</sub> = 1/√3, ρ = 2500 kg/m³).

Propagation of pressure disturbances in the host rock medium reveals that the computed pressure oscillations, originated by the ingression of buoyant magma in the magma chamber, translate into Ultra Long Period ground displacement dynamics with amplitudes of millimeter to micrometer order (Fig. 5; Longo et al. 2012b). ULP ground movements like those predicted by the present modeling could not be detected by classical broadband seismometers (although more recent seismometers extend their working range up to 100-200 s periods), while they are visible in the records from other instruments, especially
<figure><img src="image_43.png" /><figcaption>**Fig. 5** Synthetic seismic signal (blue detrended as to represent an instrumental record, black filtered [0.001,0.01] Hz) and corresponding frequency spectrum for the vertical component at all synthetic stations. Example from simulation #1</figcaption></figure>

borehole dilatometers characterized by high signal-to-noise ratio (Sacks et al. 1971).

# 4 Discussion and Conclusions

The arrival of fresh magma into an already emplaced reservoir and the consequent internal dynamics have often been invoked as possible eruption triggers, especially at Phlegraean Fields (Arienzo et al. 2010); on the other hand, footprints of magma chamber rejuvenation are often found in intrusive granites as well, testifying that it does not necessarily lead to eruption. Magma chamber rejuvenation can thus be regarded as a prototypical volcanic unrest process, that can lead to eruption or not depending on the specific system conditions (e.g. host rock compliance, volume and volatile content of injected magma).

The time scales for rejuvenation processes to be effective in magmatic reservoirs are relatively short, on the order of hours. This is consistent with what has been observed from the analyses of erupted products (Fourmentraux et al. 2012) as well as from experiments on diffusive fractionation (Perugini et al. 2010), and opens a completely new perspective in terms of unrest

duration: indications of mixing in erupted products suggest that recharge events can happen within a very short time frame from eruption, otherwise the evidence would be wiped out by the efficient mixing process (Montagna et al. 2015).

In terms of geophysical observables, convective mingling dynamics in magmatic reservoirs is associated with ultra-long-period seismic signals, characterized by frequencies in the range $10^{-2} - 10^{-4}$ Hz (Longo et al. 2012b).

The results obtained by our modeling specifically refer to the Phlegraean system. They can be extended to many other volcanoes where evidence of rejuvenation has been observed in similar magmatic settings, characterized by relatively primitive magmas that are not very different in composition and temperature, such as e.g. Mount Etna (Viccaro et al. 2006) or Stromboli (La Felice et al. 2011). For more evolved magmas, reservoirs can be dominated by crystal-rich regions (Marsh 1981; Koyaguchi and Kaneko 1999; Bachmann and Bergantz 2004; Hildreth 2004; Huber et al. 2009; Cooper and Kent 2014), and they are often at a lower temperature. The approach described above must be applied with caution in such cases, as the dynamics of the incoming primitive magma is
more likely to be described as flow through a porous medium (mush) than as fluid mingling and mixing. Nonetheless, there is some evidence for crystal-poor silicic magma reservoirs to be reactivated as well (Bachmann et al. 2002; Deering et al. 2011; Huber et al. 2012; Sliwinski et al. 2015; Wolff et al. 2015).

Given the short time scales over which the dynamical processes described here can be effective and lead to eruption, it would be beneficial to be able to routinely detect the signals described above for eruption forecasting and mitigation actions. This is especially true for long-dormant volcanoes such as Phlegraean Fields, one of the highest-risk volcanic areas in the world given the large population living within the caldera borders (Arienzo et al. 2010), for which there is still no widely accepted means of discriminating the precursors of an impending eruption (Druitt et al. 2012).

**Acknowledgements** This work has received funds from the European Union's Seventh Programme for research, technological development and demonstration under grant agreements No. 282769 VUELCO and No. 308665 MED-SUV. The manuscript has largely benefited from reviews by Fabio Arzilli and Olivier Bachmann.

## Glossary

**Magma chamber** A storage volume for magma in the crust. Typically, during their rise from the mantle magmas accumulate in regions where there are geological discontinuities.

**Rejuvenation** The process by which magmas coming from depth modify chemical and physical properties of more evolved magmas already emplaced at shallower levels.

**Primitive magma** A magmatic melt that has composition similar to that of the mantle (smaller silica content).

**Evolved magma** A magmatic melt that has undergone processes within the crust that modified its chemical properties, such as fractional crystallization and crustal assimilation.

Its composition is characterized by higher silica content.

**Convection** Exchange of mass and energy by means of cell patterns.

**Green's functions** Source to receiver transfer function; in this context through the volcanic rock medium.

## Index

- Volcanic unrest

- Magma chamber dynamics

- Magma mingling: magma mixing

- Eruption precursors

- Volcanic unrest duration

- Magma evolution

- Ground deformation

- Ultra-Long-Period seismicity

- Volcano seismicity

## References

- Aki K, Richards PG (2002) Quantitative seismology. University Science, Sausalito

- Arienzo I, Civetta L, Heumann A, Woerner G, Orsi G (2009) Isotopic evidence for open system processes within the Campanian Ignimbrite (Campi Flegrei—Italy) magma chamber. Bull Volcanol 71(3):285–300

- Arienzo I, Moretti R, Civetta L, Orsi G, Papale P (2010) The feeding system of Agnano—Monte Spina eruption (Campi Flegrei, Italy): Dragging the past into present activity and future scenarios. Chem Geol 270 (1–4):135–147

- Bachmann O, Bergantz GW (2003) Rejuvenation of the Fish Canyon magma body: a window into the evolution of large-volume silicic magma systems. Geology 31(9):789–792

- Bachmann O, Bergantz GW (2004) On the origin of crystal-poor rhyolites: Extracted from batholithic crystal mushes. J Petrol 45:1565–1582

- Bachmann O, Bergantz GW (2008) The Magma Reservoirs That Feed Supereruptions. Elements 4(1):17–21

- Bachmann O, Dungan MA, Lipman PW (2002) The Fish Canyon magma body, San Juan volcanic field, Colorado: Rejuvenation and eruption of an upper-crustal batholith. J Petrol 43:1469–1503
Bain AA, Jellinek AM, Wiebe RA (2013) Quantitative field constraints on the dynamics of silicic magma chamber rejuvenation and overturn. Contrib Min Petr 165(6):1275-1294

Cannatelli C, Lima A, Bodnar RJ, De Vivo B, Webster JD, Fedele L (2007) Geochemistry of melt inclusions from the Fondo Riccio and Minopoli 1 eruptions at Campi Flegrei (Italy). Chem Geol 237:418-432

Cooper KM, Kent AJR (2014) Rapid remobilization of magmatic crystals kept in cold storage. Nature 506:480-483

Deering CD, Bachmann O, Vogel TA (2011) The Ammonia Tanks Tuff: erupting a melt-rich rhyolite cap and its remobilized crystal cumulate. Earth Planet Sc Lett 310:518-525

De Siena L, Del Pezzo E, Bianco F (2010) Seismic attenuation imaging of Campi Flegrei: evidence of gas reservoirs, hydrothermal basins, and feeding systems. J Geophysi Res 115(B9):1-18

Di Renzo V, Arienzo I, Civetta L, D'Antonio M, Tonarini S, Di Vito MA, Orsi G (2011) The magmatic feeding system of the Campi Flegrei caldera: architecture and temporal evolution. Chem Geol 281(3-4):227-241

Druitt TH, Costa F, Deloule E, Dungan M, Scaillet B (2012) Decadal to monthly timescales of magma transfer and reservoir growth at a caldera volcano. Nature 482(7383):77-80

Elders WA, Friéleifsson GÓ, Zierenberg RA, Pope EC, Mortensen AK, Guðmundsson Á, Lowenstern JB, Marks NE, Owens L, Bird DK, Reed M (2011) Origin of a rhyolite that intruded a geothermal well while drilling at the Krafla volcano Iceland. Geology 39(3):231-234

Fourmentraux C, Metrich N, Bertagnini A, Rosi M (2012) Crystal fractionation, magma step ascent, and syn-eruptive mingling: the Averno 2 eruption (Phlegrean Fields, Italy). Contrib Min Petr 163(6):1121-1137

Giordano D, Russell JK, Dingwell DB (2008) Viscosity of magmatic liquids: a model. Earth Planet Sci Lett 271:123-143

Girard G, Stix J (2009) Magma recharge and crystal mush rejuvenation associated with early post-collapse Upper Basin Member rhyolites, Yellowstone caldera. Wyoming. J Petrol 50(11):2095-2125

Gudmundsson A (2006) How local stresses control magma-chamber ruptures, dyke injections, and eruptions in composite volcanoes. Earth Sci Rev 79:1-31

Hildreth WS (2004) Volcanological perspectives on Long Valley, Mammoth Mountain, and Mono Craters: several contiguous but discrete systems. J Vol Geoth Res 136:169-198

Huber C, Bachmann O, Manga M (2009) Homogenization processes in silicic magma chambers by stirring and mushification (latent heat buffering). Earth Planet Sci Lett 283:38-47

Huber C, Bachmann O, Dufek J (2012) Crystal-poor versus crystal-rich ignimbrites: a competition between stirring and reactivation. Geology 40:115-118

Koyaguchi T, Kaneko K (1999) A two-stage thermal evolution model of magmas in continental crust. J Petrol 40:241-254

La Felice S, Landi P (2011) The 2009 paroxysmal explosions at Stromboli (Italy): magma mixing and eruption dynamics. Bull Volc 73(9):1147-1154

Longo A, Barsanti M, Cassioli A, Papale P (2012a) A finite element Galerkin/least-squares method for computation of multicomponent compressible in compressible flows. Comp Fluids 67:57-71

Longo A, Papale P, Vassalli M, Saccorotti G, Montagna CP, Cassioli A, Giudice S, Boschi E (2012b) Magma convection and mixing dynamics as a source of Ultra-Long-Period oscillations. Bull Volcanol 74:873-880

Mangiacapra A, Moretti R, Rutherford MJ, Civetta L, Orsi G, Papale P (2008) The deep magmatic system of the Campi Flegrei caldera (Italy). Geophys Res Lett 35

Marsh BD (1981) On the crystallinity, probability of occurrence, and rheology of lava and magma. Contrib Min Petr 78:85-98

Mollo S, Masotta M, Forni F, Bachmann O, De Astis G, Moore G, Scarlato P (2015) A K-feldspar-liquid hygrometer specific to alkaline differentiated magmas. Chem Geol 392:1-8

Montagna CP, Papale P, Longo A (2015) Timescales of mingling in shallow magmatic reservoir. In: Caricchi L, Blundy J D (eds) Chemical, physical and temporal evolution of volcanic systems. Geological Society, London, Special Publications 422

O'Brien GS, Bean CJ (2004) A 3D discrete elastic lattice method for seismic wave propagation in heterogeneous media with topography. Geophys Res Lett 31, L14608

Papale P, Moretti R, Barbato D (2006) The compositional dependence of the multicomponent volatile saturation surface in silicate melts. Chem Geol 229:78-95

Pappalardo L, Ottolini L, Mastrolorenzo G (2007) The Campanian Ignimbrite (southern Italy) geochemical zoning: insight on the generation of a super-eruption from catastrophic differentiation and fast withdrawal. Contrib Min Petr 156:1-26

Perugini D, Poli G, Petrelli M, Campos CP, Dingwell DB (2010) Time-scales of recent Phlegrean Fields eruptions inferred from the application of a diffusive fractionation model of trace elements. Bull Volcanol 72(4):431-447

Petrelli M, Perugini D, Poli G (2011) Transition to chaos and implications for time-scales of magma hybridization during mixing processes in magma chambers. Lithos 125(1-2):211-220

Sacks IS, Selwyn S, Evertson DW (1971) Sacks-Evertson strainmeter, its installation in Japan and some preliminary results concerning strain steps. P Jpn Acad 47(9):707-712
Signorelli S, Vaggelli G, Carroll C, Romano MR (2001) Volatile element zonation in Campanian Ignimbrite magmas (Phlegrean Fields, Italy): evidence from the study of glass inclusions and matrix glasses. Contrib Min Petr 140:543–553

Sliwinski JT, Bachmann O, Ellis BS, Dávila-Harris P, Nelson BK, Dufek J (2015) Eruption of shallow crystal cumulates during caldera-forming events on Tenerife. Canary Islands, J Petrol 56(11):2173–2194

Sparks R, Sigurdsson H, Wilson L (1977) Magma mixing: a mechanism for triggering acid explosive eruptions. Nature 267:315–318

Till CB, Vazquez JA, Boyce JW (2015) Months between rejuvenation and volcanic eruption at Yellowstone caldera. Wyoming. Geology 43(8):695–698

Viccaro M, Ferlito C, Cortesogno L, Cristofolini R, Gaggero L (2006) Magma mixing during the 2001

event at Mount Etna (Italy): Effects on the eruptive dynamics. J Vol Geo Res 149(1–2):139–159

Voight B, Widiwijayanti C, Mattioli G S, Elsworth D, Hidayat D, Strutt M (2010) Magma-sponge hypothesis and stratovolcanoes: Case for a compressible reservoir and quasi-steady deep influx at Soufriere Hills Volcano, Montserrat. Geophys Res Lett 37(19): L00E05

Wolff JA, Ellis BS, Ramos FC, Starkel WA, Boroughs S, Olin PH, Bachmann O (2015) Remelting of cumulates as a process for producing chemical zoning in silicic tuffs: A comparison of cool, wet and hot, dry rhyolitic magma systems. Lithos 236–237:275–286

Zollo A, Maercklin N, Vassallo M, Dello Iacono D, Virieux J, Gasparini P (2008) Seismic reflections reveal a massive melt layer feeding Campi Flegrei caldera. Geophys Res Lett 35(12):L12306

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.
# Magma Mixing: History and Dynamics of an Eruption Trigger

Daniele Morgavi, Ilenia Arienzo, Chiara Montagna, Diego Perugini and Donald B. Dingwell

## Abstract

The most violent and catastrophic volcanic eruptions on Earth have been triggered by the refilling of a felsic volcanic magma chamber by a hotter more mafic magma. Examples include Vesuvius 79 AD, Krakatau 1883, Pinatubo 1991, and Eyjafjallajökull 2010. Since the first hypothesis, plenty of evidence of magma mixing processes, in all tectonic environments, has accumulated in the literature allowing this natural process to be defined as fundamental petrological processes playing a role in triggering volcanic eruptions, and in the generation of the compositional variability of igneous rocks. Combined with petrographic, mineral chemistry and geochemical investigations, isotopic analyses on volcanic rocks have revealed compositional variations at different length scales pointing to a complex interplay of fractional crystallization, mixing/mingling and crustal contamination during the evolution of several magmatic feeding systems. But to fully understand the dynamics of mixing and mingling

<footnote>D. Morgavi · D. Perugini
Department of Physics and Geology,
University of Perugia, Piazza Università,
06100 Perugia, Italy</footnote>

<footnote>D. Morgavi (✉) · D.B. Dingwell
Department Earth and Environmental Sciences,
Ludwig-Maximilian-University (LMU),
Theresienstrasse 41/III, 80333 Munich, Germany
e-mail: daniele.morgavi@unipg.it</footnote>

<footnote>I. Arienzo
Istituto Nazionale Di Geofisica E Vulcanologia,
Osservatorio Vesuviano, Via Diocleziano, 328,
8124 Naples, Italy</footnote>

<footnote>C. Montagna
Istituto Nazionale Di Geofisica E Vulcanologia,
Via Della Faggiola, 32, Pisa 56126, Italy</footnote>
processes, that are impossible to observe directly, at a realistically large scale, it is necessary to resort to numerical simulations of the complex interaction dynamics between chemically different magmas.

Keywords

Magma mixing • Mingling • Isotope • Modelling

# 1 Magma Mixing: A Brief Historical Overview

One of the first investigations on magma mixing recorded in the literature is the work of the chemist Bunsen (1851), a scholar at the University of Heidelberg who published research on the chemical variation of igneous rocks from the western region of Iceland. Through chemical analyses Bunsen highlighted that the linear correlation between pairs of chemical elements in binary plots in those Icelandic rocks was the consequence of “simple” binary mixing between two magmas with different chemical composition. Bunsen published this data and, for the first time, magma mixing was taken into account to explain the chemical variation of a suite of igneous rocks. This idea triggered a strong critical reaction from the geological community; the strongest opposition coming from Wolfgang Sartorius Freiherr von Waltershausen an expert on the Iceland and Etna volcanic areas at that time. He mostly argued against the method used by Bunsen of averaging rock analyses to calculate the starting end-members that eventually took part in the mixing process. Sartorius criticized not only the arbitrary choice of the end-members but also disliked the idea of Bunsen of an extensive layering of felsic/mafic rocks and magmas beneath Iceland.

Since the beginning of the 20th century, the experimental and thermodynamic work of Norman L. Bowen (e.g., Bowen 1928) has had a profound influence upon the way petrological

processes and igneous differentiations are conceived. The conceptual model of fractional crystallization was firmly established as the most fundamental petrological process for generating the diversity of igneous rocks and remained so for many decades. Although Bowen did not explicitly deny the possibility of magma mixing, he reinterpreted field evidence of magma mixing rather as immiscibility of liquids (e.g., Bowen 1928). In 1944, Wilcox published a work on the Gardner River complex (Yellowstone, USA; Wilcox 1944) which is now considered a milestone for evidence of magma mixing, even if at that time it received strong comments from Fenner and remained one of the few papers on the topic. Only in the 1970's geoscientists started to deeply investigate magma mixing, recorded as a plethora of unequivocal evidence in both plutonic and volcanic rocks, as a major petrogenetic process (e.g., Eichelberger 1978, 1980; Blundy and Sparks 1992; Wiebe 1994; Wilcox 1999). Since the first hypotheses about the origin of mixed igneous rocks (e.g., Bunsen 1851), plenty of evidence of magma mixing processes, in all tectonic environments, throughout geological time, has accumulated in the literature allowing this natural process to be defined as a fundamental petrological process playing a key role in the generation of the compositional variability of igneous rocks and as a major process for planetary differentiation (e.g., Eichelberger 1978, 1980; Blundy and Sparks 1992; Wiebe 1994; De Campos et al. 2004; Perugini and Poli 2012; Morgavi et al. 2016).
## 2 Magma Mixing: Field Evidence

It is common practice in the petrological community to split magma mixing into two separate physico/chemical processes: (i) mechanical mixing (also referred to as “magma mingling”), by which two or more batches of magma mingle without chemical exchanges between them, and (ii) a chemical mixing (also referred to as “magma mixing”) triggered by chemical exchanges between the interacting magmas in which elements move from one magma to the other according to compositional gradients continuously generated by the mechanical dispersion of the two magmas (e.g., Flinders and Clemens 1996). Physically, “magma mingling” is mainly controlled by the viscosity contrast between the two magmas; decreasing of the viscosity contrast results in progressively more efficient mingling dynamics (e.g., Sparks and Marshall 1986; Grasset and Albarede 1994; Bateman 1995; Poli et al. 1996; Perugini and Poli 2005). Chemically, “magma mixing” is driven by the mobility of chemical elements in the melt fractions of the two magmas (e.g., Lesher 1990; Baker 1990). Linear variations in inter-elemental plots for a set of rock samples have long been considered as the sole evidence for the occurrence of magma mixing (e.g., Fourcade and Allegre 1981).

The adoption of the above conceptual models led to the common practice of applying the term magma mingling to indicate the process acting to physically disperse (no chemical exchanges are involved) two or more magmas, whereas the term magma mixing indicates that the mingling process is also accompanied by chemical exchanges. Although such a conceptual approach may allow us to simplify the complexity of the magma mixing process and make it more tractable from the petrological point of view, unfortunately such terminology is not consistently used in the literature and this causes some misunderstanding. Although it is not always easy to clearly discriminate between the two processes, mingling is quite a rare process in nature as physical dispersion and chemical exchanges must occur in tandem during magma mixing processes (e.g., Wilcox 1999; Perugini and Poli 2012).

The most common evidence for magma mixing in igneous rocks is the occurrence of textural heterogeneity; the processes responsible for this have been discussed extensively in many works in the last decades (e.g., Eichelberger 1975; Anderson 1976; Bacon 1986; Didier and Barbarin 1991; Wada 1995; De Rosa et al. 1996; Ventura 1998; Smith 2000; Snyder 2000; De Rosa et al. 2002; Perugini et al. 2002, 2007; Perugini and Poli 2005, 2012; Pritchard et al. 2013; Morgavi et al. 2016).

In order to provide possible classification of magma mixing structures, the evidence of mechanical mixing in igneous rocks can be roughly divided into three different groups: (i) flow structures, (ii) magmatic enclaves and (iii) physico-chemical disequilibria in melts and crystals (e.g., Walker and Skelhorn 1966; Didier and Barbarin 1991; Hibbard 1995; Flinders and Clemens 1996; Wilcox 1944, 1999; Perugini et al. 2002, 2003; Streck 2008; Perugini and Poli 2012; Morgavi et al. 2016). Flow structures can be readily recognized in field outcrops as they show alternating light and dark coloured bands constituted by magmas with different compositions. Figure 1a, b shows some examples of fluid structures occurring in volcanic rocks from Grizzly Lake outcrop in Yellowstone National Park (USA) (Pritchard et al. 2013; Morgavi et al. 2016) and from Soufrière Hills volcano (Island of Montserrat, UK) (Plail et al. 2014). In particular, Fig. 1a shows flow bands of rhyolitic magma (white) intruding in a basaltic/hybrid magma (red to dark grey) whereas Fig. 1b shows an alternation of flow bands of rhyolitic (white) and hybrid magma (dark grey) across which basaltic enclaves (light grey) occur. The latter are surrounded by flow bands of hybrid filaments (blue).

Magmatic enclaves are probably the structural evidence that, according to common thinking, mostly characterize magma mixing processes. The term magmatic enclave is used to identify a discrete portion of a magma occurring within a host magma with a different composition (e.g., Wilcox 1944; Walker and Skelhorn 1966; Bacon 1986, Didier and Barbarin 1991). Generally, enclaves display quite sharp contacts with the
<figure><img src="image_44.png" /><figcaption>**Fig. 1** Detailed images of the mixing features present at Grizzly Lake (Yellowstone) (a, b) and Soufrière Hills volcano (Montserrat) (c, d). Figure 1a shows the rhyolitic magma (white) has apparently intruded into the basaltic magma and the hybrid portions are present at the contact between the two end-members. Two large basaltic enclaves are visible at the bottom left and at the centre</figcaption></figure>

right. Figure 1b shows the stretching and folding of hybrid magma (dark grey) into a rhyolitic portion (white) with the presence of several basaltic enclaves. Figure 1c shows a basaltic enclave surrounded by andesitic magma (Soufrière Hills, from the 2010 eruption). Figure 1d exhibits at the centre a basaltic enclave in and andesitic host from Soufrière Hills volcano

host rock, although it is not rare to observe that some enclaves display engulfment and disruption of their boundaries due to infiltration of the host magma. Some examples of enclaves found in the volcanic rocks from Soufrière Hills are shown in Fig. 1c, d.

Disequilibrium textures in minerals (Fig. 2a–c) can be viewed as recorders of the thermal and compositional disequilibria operating in the magmatic system during the development of magma mixing processes. As the zoning pattern can be well preserved in minerals from both the plutonic and volcanic rocks, crystal populations from both environments can be used

to reconstruct the time evolution of thermal and compositional exchanges between the two magmas during mixing. Recent studies highlighted the importance of detailed investigations of crystal compositional variability not only to reconstruct the fluid-dynamic regime governing the evolution of the igneous body, but also to understand the length-scale of the compositional variability induced by the mixing process, the latter being considered as a proxy to estimate the residence time of magmas in sub-volcanic reservoirs prior to eruption (e.g., Costa and Chakraborty 2004; Martin et al. 2008; Chamberlain et al. 2014; Perugini et al. 2003).
<figure><img src="image_45.png" /></figure>

**Fig. 2** Backscattered electron (BSE) image of a section of disequilibrium textures in rock and minerals from Yellowstone, Montserrat and Phlegrean Fields. a Mixed rock from Grizzly Lake Complex (Yellowstone) showing the interaction between the basaltic portion, the hybrid portion and the rhyolitic portion. b Mixed rock section from the 2010 Soufrière Hill eruption (Montserrat) showing disequilibrium texture in the basaltic and the

andesitic rock. c Crystal from the 4.67 ± 0.09 cal ka Agnano Monte-Spina eruption (Phlegrean fields) occurred from a vent in the Agnano-San Vito area, has a darker (i.e., Ba-poorer) resorbed core and an inner rim with “swirly” zonation textures that indicate crystallization and dissolution. The outer rim is characterized by small-scale wavy oscillatory zoning that results from high-frequency growth and resorption events

## 3 Numerical and Experimental Studies: New Ideas for Deciphering the Complexity of Magma Mixing

Studies focused on numerical and experimental investigation of magma mixing dynamics (e.g., Perugini et al. 2003, 2008, 2015; De Campos et al. 2004, 2008, 2011; Petrelli et al. 2011; Montagna et al. 2015; Morgavi et al. 2015;

Laumonier et al. 2014) can provide additional tools for a better understanding of the complexity of the mixing process, the evolution of which is governed by a continuous exchanges. One of the most striking results arising from these studies is that, during mixing, chemical elements experience a diffusive fractionation process due to the development in time of chaotic mixing dynamics (Perugini et al. 2006, 2008). This process is considered the source of strong deviations in many chemical elements from the linear
variations in inter-elemental plots that would otherwise be expected, based on a conceptual model classically adopted in the geochemical modelling of magma mixing processes (e.g., Fourcade and Allegre 1981; Perugini and Poli 2012 and references therein). Recent studies on the mineralogical and geochemical features of mixed rocks (e.g., Hibbard 1981, 1995; Wallace and Bergantz 2002; Costa and Chakraborty 2004; Perugini and Poli 2005; Slaby et al. 2010), as well as those focused on quantitative analyses of morphologies related to textural heterogeneity (e.g., Wada 1995; De Rosa et al. 2002; Perugini and Poli 2005; Perugini et al. 2002, 2003) have highlighted the dominant role played by chaotic mixing dynamics in producing the substantial complexity of geochemical variations and textural patterns found in the resultant rocks (e.g., Flinders and Clemens 1996; De Campos et al. 2011; Morgavi et al. 2013a, b, c, 2016). Despite significant attention in the past, however, few works have focused on the understanding of the relationship between the morphology of the mixing patterns and the geochemical variability of the system using experimental devices (e.g., De Rosa et al. 2002).

Based upon the combination of field observations and the outcome from numerical simulations, a new experimental apparatus has been developed to perform mixing experiments using high viscosity silicate melts at high temperature (De Campos et al. 2011; Morgavi et al. 2013a, c, 2015). This device has been used to study the mixing process between natural melts, enabling the investigation of the influence of chaotic dynamics on the geochemical evolution of the system of mixing magmas (Morgavi et al. 2013a, b, c, 2015).

Preliminary results indicate that the time evolution of compositional exchanges between magmas from the experiments can be effectively modelled, leading to the prospect that the record of magma mixing processes may serve as chronometers to estimate the time interval between mixing and eruption (Perugini et al. 2010; Perugini and Poli 2012; Morgavi et al. 2013a, b, c, 2015; Perugini et al. 2015).

## 4 Geochemical Evidence of Magma Mixing/Mingling: An Example from the Campi Flegrei Volcanic Area

In some volcanic areas chemically and isotopically distinct magmas have been erupted, and their composition identified by analyzing the chemical composition of the erupted products (e.g., Pantelleria (Italy), Gedemsa and Fanta 'Ale (Main Ethiopian Rift), Gorely Eruptive Center (Kamchatka); Civetta et al. 1997; Giordano et al. 2014; Seligman et al. 2014). However, in other volcanic complexes the majority of the erupted products are chemically rather homogeneous, displaying a dominant composition (e.g., trachybasalt at Mt. Etna; trachyte at Campi Flegrei, Italy; basalt at Réunion Island, Indian Ocean; andesite at Popocatepetl, Mexico). Despite the roughly homogenous composition of products from these volcanic areas, their isotopic features suggest that complex open system processes occurred and superposed the main fractional crystallization trend. In fact, isotopic analyses (e.g., Sr, Nd, Pd, B) have been proven to be an important tool for discriminating between closed-system fractional crystallization and open-system magma mixing/mingling or crustal contamination (e.g., James 1982; Knesel et al. 1999; Turner and Foden 2001). Combined with petrographic, mineral chemistry and chemical investigations, isotopic analyses on volcanic rocks have revealed compositional variations at different length-scales (bulk rock, minerals, single crystals) pointing to a complex interplay of fractional crystallization, mixing/mingling and crustal contamination during the evolution of several magmatic feeding systems (e.g., Di Renzo et al. 2011 and references therein; Melluso et al. 2012; Corsaro et al. 2013 and references therein; Di Muro et al. 2014 and references therein; Brown et al. 2014). Furthermore, together with conventional isotopic analyses, current technologies permit high precision, in situ determination of Sr isotopic ratios of portions of phenocryst and glasses. In fact, microsampling by MicroMill™, coupled with isotopic measurement by Thermal
Ionization Mass Spectrometry (TIMS), allows for the performance of high precision determination of Sr isotopic composition of single crystals or portions of them. This information, unobtainable from bulk samples, has been used successfully to gather information on the time- and length-scales of the pre-eruptive magmatic processes, for identifying mantle sources and/or magmatic end-members and for tracking the time evolution of magma differentiation (e.g., Davidson et al. 1990; Davidson and Tepley 1997; Davidson et al. 1998; Knesel et al. 1999; Font et al. 2008; Kinman et al. 2009; Francalanci et al. 2012; Braschi et al. 2012; Jolis et al. 2013; Arienzo et al. 2015).

Among the active volcanic areas worldwide the volcanic hazard posed by the Campi Flegrei caldera is extremely high, due to its explosive character. Both the high volcanic hazard and the intense urbanization result in an extreme volcanic risk in this area, leading to a considerable interest in understanding which processes might contribute to triggering of eruptions and controlling? eruptive dynamics. The Campi Flegrei caldera is a nested and resurgent structure in the Campania Region, South Italy (Orsi et al. 1996), possibly formed after two large caldera forming eruptions: the Campanian Ignimbrite eruption (39 ka, Fedele et al. 2008) and the Neapolitan Yellow Tuff (15 ka, Deino et al. 2004). Its magmatic system is still active as testified by the occurrence of the last eruption in 1538 AD, as well as the present widespread fumaroles and hot springs activity, and the persistent state of unrest (Del Gaudio et al. 2010; Chiodini et al. 2003, 2012, 2015; Moretti et al. 2013). For compositionally homogenous magmas such as those extruded at the Campi Flegrei caldera (trachytes and phonolites being by far the most abundant rocks), major oxide and trace element variations cannot be used to unequivocally establish which magma evolution processes operated. Thus, together with petrographic, mineral chemistry and chemical data, isotopic investigations on volcanic rocks spanning the history of the volcano have been performed in recent decades in order to define the role of variable magmatic processes in

the evolution of its feeding system up to eruption (e.g., Civetta et al. 1997; D'Antonio et al. 1999, 2013; de Vita et al. 1999; Pappalardo et al. 2002; Fedele et al. 2008, 2009; Tonarini et al. 2004, 2009; Arienzo et al. 2010, 2011; Perugini et al. 2010; Di Vito et al. 2011; Melluso et al. 2012; Arienzo et al. 2015).

In particular, detailed investigations of the geochemical and isotopic (Sr, Nd, Pb, and B) features of the younger than 15 ka Campi Flegrei volcanic products gave understanding to how many variable magmatic components, rising from large depth and/or stagnating in middle crustal reservoir(s), recharged the shallowest reservoir(s) and interacted with magma batches left from previous eruptions (Di Renzo et al. 2011). One identified magmatic component, geochemically similar to magma from the Neapolitan ca. 0.70750, 143Nd/144Nd ratio of ca. 0.51247, 206Pb/204Pb of ca. 19.04 and d11B of ca. -7.8‰), has been the most prevalent component over the past 15 ka. A second magmatic component, having geochemical features similar to the Minopoli 2 magma (D'Antonio et al. 1999; Di Renzo et al. 2011), first erupted 10 ka ago, is shoshonitic in composition. It is the most enriched in radiogenic Sr (87Sr/86Sr of ca. 0.70850) and unradiogenic Nd and Pb (143Nd/144Nd ratio of ca. 0.51238, 206Pb/204Pb of ca. 18.90), and it is characterised by the lowest d11B value of ca. -7.4‰. The third component is trachytic in composition and is characterized by lower 206Pb/204Pb (ca. 19.08), 87Sr/86Sr (ca. 0.70720) and d11B (-9.8‰) and higher 143Nd/144Nd (ca. 0.51250), with respect to the Neapolitan Yellow Tuff component (Tonarini et al. 2009; Di Renzo et al. 2011; Arienzo et al. 2015). This third composition is known as the Astroni 6 component due to the fact that it best recognized in the Astroni 6 erupted products (Di Renzo et al. 2011). During the past 5 ka of activity, this new component has been suggested to have mixed in variable proportions with the Neapolitan Yellow Tuff and Minopoli 2 magmatic components, which dominated the Campi Flegrei volcanic activity mostly in the time span from 15 to 5 ka (Fig. 3).
<figure><img src="image_46.png" /></figure>

**Fig. 3** Modelling of the Sr-Nd isotopic features of some of the Campi Flegrei volcanics of the past 5 ka, by assuming mixing among the Astroni 6 (Ast-6)-, Neapolitan Yellow Tuff (NYT)- and Minopoli 2 (Min 2)-like magmatic components. The green, yellow and black boxes represent the range of Sr and Nd isotopes of the products erupted during the Astroni 6, Neapolitan Yellow

Tuff and the Minopoli 2 eruptions, respectively. Symbols inside the plot represent volcanic products belonging to the listed eruptions. The vertical error bar is the uncertainty in 143Nd/144Nd determination at the 2σ level of confidence; that for 87Sr/86Sr is included in the symbols. Modified after Arienzo et al. (2016)

Based on isotope investigations and melt inclusions studies, Arienzo et al. (2016) suggested that the Astroni 6 component, although undergoing differentiation during uprising, had a deep origin (larger than 8 km depth). Indeed, this magma rose not only inside and along the margins of the caldera, but also at the intersection between SE-NW and NE-SW regional fault systems mixing with the NYT-like magma component at shallower depth, and possibly entrapping crystals accumulated during older eruptions.

This detailed study of the Campi Flegrei volcanic system highlights that Sr isotopic microanalysis and, in general, more conventional isotopic analyses, coupled with petrographic, mineral chemistry and geochemical data can provide a better knowledge of the mixing/mingling processes and of the mixing end-members. In turn, they provide (i) information for evaluating the

volcanic hazards and mitigating the related risks and (ii) the basic geochemical and petrologic knowledge for the numerical simulations.

**5 Numerical Simulation of Magma Mingling and Mixing**

To understand mixing and mingling processes at a realistically large scale, it is necessary to resort to numerical simulations of the complex interaction dynamics between chemically different magmas. Referring to the archetypal case of the Campi Flegrei magmatic system as described above, the interaction of a shoshonite and a more evolved phonolite has been investigated in detail to provide constraints on the time and length scales of the mixing dynamics. The simulated system consists of a very large, deep (8 km) reservoir connected by a dike to a shallower,
<figure><img src="image_47.png" /><figcaption>**Fig. 4** Initial setup for the numerical simulations of magma chamber replenishment at Campi Flegrei caldera. The two interacting magmas are characterized by different compositions and volatile contents, thus densities and gas volume fractions</figcaption></figure>

smaller chamber (Fig. 4). The chemical interactions between the two magmas cannot be resolved on the simulated large scale, as the computational costs required would be too high. The shoshonitic magma, being richer in volatiles than the resident phonolitic melt, rises into the phonolitic chamber by buoyancy, generating convection and mixing within the reservoir and in the feeding dike. Both magmatic components include a liquid (silicate melt) and a gaseous phase, that cannot decouple from the host melt. Space-time varying volatile exsolution is computed as a function of local composition, pressure and temperature following Papale et al. (2006). More details on magma chamber dynamics can be found in Chap. 8 ‘Magma chamber rejuvenation: insights from numerical models’ of this book.

The main results show that the chaotic patterns observed in the products and in recently developed experimental setups (Morgavi et al. 2013a, b, c, 2015) are reproducible (Fig. 5), and

that the two magmas mingle very efficiently from the beginning of their interaction. As time progresses, convection slows down due to smaller buoyancy of the incoming mixed component, and the instability proceeds in time asymptotically: the more the two end-members have mingled, the less intense the convection. A time-dependent mixing efficiency $\eta_C$ can be defined as:

$$ \eta_C = \frac{|m_R(t) - m_R(0)|}{m_R(0)}. \quad (1) $$

In the equation above, $m_R(t)$ is the mass of the resident phonolitic magma at time $t$, thus $m_R(0)$ is the initial phonolite mass. The mixing efficiency $\eta_C$ represents the relative variation of the mass of the initially resident magma in a certain region of the domain. Figure 6 shows the time evolution of mixing efficiency in the shallow chamber, for different simulated setups in terms of chamber geometry and volatile content. It clearly shows
<figure><img src="image_48.png" /><figcaption>**Fig. 5** Magma mixing in natural samples, experimental setups 684 and numerical modelling. Left: Lesvos (Greece) lava flow, from Perugini and Poli (2012); centre:</figcaption></figure>

experimental setup, from Morgavi et al. (2013b); right: simulations of magma chamber replenishment

<figure><img src="image_49.png" /><figcaption>**Fig. 6** Variation of convection efficiency with time in the shallow chamber for different simulated scenarios, characterized by varying geometry and total volatile content in the shallow chamber</figcaption></figure>

that mingling is very effective on relatively short time scales, on the order of hours, in agreement with mixing timescales derived from the geochemical modelling of magma mixing experiments (Perugini et al. 2010). When buoyancy drives fast dynamics, more than 40% of the original end-member magmas have mingled in the feeder dyke within 4 h from the arrival of the gas-rich magma from depth. The asymptotic behaviour seems to be reached earlier for less efficient setups: after 4–5 h from the onset of the instability, the system seems to have reached a quasi-steady state. In these cases, a much smaller part of the two magmas has mingled.

## 6 Magma Mixing Time Scale and Eruption Trigger

The most violent and catastrophic volcanic eruptions on Earth have been triggered by the refilling of a felsic volcanic magma chamber by a hottest and more mafic magma (Kent et al. 2010; Murphy 1998). Examples include Vesuvius 79 AD (Cioni et al. 1995), Krakatau 1883 (Self 1992), Pinatubo 1991 (Kress 1997), the Canadian Ignimbrite (Arienzo et al. 2009) and Eyjafjallajökull 2010 (Sigmundsson et al. 2010). Injection of the more mafic magma into the felsic
magma triggers convection dynamics and wide-
spread mixing (Sparks et al. 1977). Vesiculation
induced by convection increases magma pressure
and may fracture the volcanic edifice triggering
an explosive eruption. The injection and mixing
process is accompanied by geophysical signals,
such as earthquakes, gravity changes and ultra-
long-period ground oscillations, that can now be
accurately detected (Williams and Rymer 2002;
Longo et al. 2012; Bagagli et al. 2017). The
knowledge of the time elapsing between the
beginning of mixing (and associated geophysical
signals) and eruption is thus of greatest impor-
tance in forecasting the onset of a volcanic
eruption.

Recent studies highlighted that in order to
preserve magma mixing structures (i.e., fila-
ments, swirls, bandings) in the rocks, the time
elapsed between the beginning of mixing and the
subsequent eruption must be very short; on the
order of hours or days (Perugini et al. 2010).
Preserved structures would indicate mixing
was the last process recorded by the magmatic
system and its study can unravel unprecedented
information on pre-eruptive behaviour of
volcanoes.

The compositional heterogeneity produced by
magma mixing, and subsequently frozen in time
in the volcanic products, can hence be viewed as
a broken clock at a crime scene; it can potentially
be used to determine the time of the incident.
Following this idea and combining numerical
simulations with magma mixing experiments
using natural compositions and statistical analy-
ses it was shown that for three volcanic eruptions
from the Campi Flegrei volcanic system
(Astroni, Averno and Agnano Monte Spina) the
mixing-to-eruption timescale are of the order of a
few minutes (Perugini et al. 2015). These time-
scales indicate that very little time elapsed from
the moment mixing started until eruption. These
results are in agreement with recent numerical
simulations of magma mixing (Montagna et al.
2015) that highlight mixing timescales of a few
hours to attain complete hybridization of mag-
mas for the Campi Flegrei magmatic systems.

These results have implications for civil pro-
tection planning of future volcanic crisis as the

high velocities of ascending magmas may imply
little warning time in volcanic crises. These
findings can be a starting point towards a uni-
fying model explaining chemical exchanges in
magmatic systems and supplying information on
the use of chemical element mobility as
geochronometers for volcanic eruptions. This
may provide unparalleled clues for building an
inventory of past and recent volcanic eruption
timescales and could be decisive for hazard
assessment in active volcanic areas.

**Acknowledgements** This research was funded by the European Union's Seventh Programme for research technological development and demonstration under grant agreement No 282759—VUELCO and by the ERC Consolidator Grant 612776—CHRONOS.

References

Anderson AT (1976) Magma mixing: petrological process and volcanological tool. J Volcanol Geotherm Res 1:3–33

Arienzo I, D'Antonio M, Di Renzo V, Tonarini S, Minolfi G, Orsi G, Carandente A, Belviso P, Civetta L (2015) Isotopic microanalysis sheds light on the magmatic endmembers feeding volcanic eruptions: the Astroni 6 case study (Campi Flegrei, Italy) J Volcanol Geoth Res 304:24–37

Arienzo I, Mazzeo FC, Moretti R, Cavallo A, D'Antonio M (2016) Open-system magma evolution and fluid transfer at Campi Flegrei caldera (Southern Italy) during the past 5 ka as revealed by geochemical and isotopic data: the prototype of Nisida eruption. Chem Geol 427:109–124

Bacon CR (1986) Magmatic inclusions in silicic and intermediate volcanic rocks. J Geophys Res 91:6091–6112

Bagagli M, Montagna CP, Papale P, Longo A (2017) Signature of magmatic processes in strainmeter records at Campi Flegrei (Italy). Geophys Res Lett 44(2):718–725

Baker DR (1990) Chemical Interdiffusion of Dacite and Rhyolite—Anhydrous Measurements at 1 Atm and 10 Kbar, Application of Transition-State Theory, and Diffusion in Zoned Magma Chambers. Contrib Mineral Petr 104:407–423

Bateman R (1995) The interplay between crystallization, replenishment and hybridisation in large felsic magma chambers. Earth Sci Rev 39:91–106

Blundy J, Sparks SRJ (1992) Petrogenesis of mafic enclaves in granitoids of the Adamello massif, Italy. J Petrol 33:1039–1104

Bowen N (1928) The evolution of the igneous rocks. Princeton University Press, Princeton, N.J
Braschi E, Francalanci L, Vougioukalakis GE (2012) Inverse differentiation pathway by multiple mafic magma refilling in the last magmatic activity of Nisyros Volcano, Greece. Bull Volcanol 74:1083–1100

Brown RJ, Civetta L, Arienzo I, D'Antonio M, Moretti R, Orsi G, Tomlinson EL, Albert PG, Menzies MA (2014) Geochemical and isotopic insights into the assembly, evolution and disruption of a magmatic plumbing system before and after a cataclysmic caldera Dcollapse eruption at Ischia volcano (Italy). Contrib Mineral Petrol 168, 1035. doi:https://doi.org/10.1007/s00410-014-1035-1

Bunsen W (1851) Über die Processe der vulkanischen Gesteinsbildungen Islands. Ann Phys Chem (Dritte Reihe) 83:197–272

Chamberlain KJ, Morgan DJ, Wilson CJN (2014) Timescales of mixing and mobilisation in the Bishop Tuff magma body: perspectives from diffusion chronometry. Contrib Mineral Petr 167:1034. https://doi.org/10.1007/s00410-014-1034-2

Chiodini G, Caliro S, De Martino P, Avino R, Gherardi F (2012) Early signals of new volcanic unrest at Campi Flegrei caldera? Insights from geochemical data and physical simulations. Geology 40:943–946

Chiodini G, Todesco M, Caliro S, Del Gaudio C, Macedonio G, Russo M (2003) Magma degassing as a trigger of bradyseismic events; the case of Phlegrean Fields (Italy). Geophys Res Lett 30(8):1434. https://doi.org/10.1029/2002GL01679

Chiodini G, Vandemeulebrouck J, Caliro S, D'Auria L, De Martino P, Mangiacapra A, Petrillo Z (2015) Evidence of thermal-driven processes triggering the 2005–2014 unrest at Campi Flegrei caldera. Earth Planet Sci Lett 414:58–67

Cioni R, Civetta L, Marianelli P, Metrich N, Santacroce R, Sbrana A (1995) Compositional layering and syn-eruptive mixing of a periodically refilled shallow magma chamber: the AD 79 Plinian eruption of Vesuvius. J Petrol 36:739–776

Civetta L, Orsi G, Pappalardo L, Fisher RV, Heiken G, Ort M (1997) Geochemical zoning, mingling, eruptive dynamics and depositional processes—the Campanian Ignimbrite, Campi Flegrei Caldera, Italy. J Volcanol Geotherm Res 75:183–219

Corsaro RA, Di Renzo V, Di Stefano S, Miraglia L, Civetta L (2013) Relationship between magmatic processes in the plumbing system of Mt. Etna and the dynamics of the eastern flank: inferences from the petrologic study of the products erupted from 1995 to 2005. J Volcanol Geotherm Res 251:75–89

Costa F, Chakraborty S (2004) Decadal time gaps between mafic intrusion and silicic eruption obtained by chemical zoning patterns in olivine. Earth Planet Sci Lett 227:517–530

D'Antonio M, Civetta L, Orsi G, Pappalardo L, Piochi M, Carandente A, de Vita S, Di Vito MA, Isaia R (1999) The present state of the magmatic system of the Campi Flegrei caldera based on a reconstruction of its behavior in the past 12 ka. J Volcanol Geotherm Res 91:247–268

D'Antonio M, Tonarini S, Arienzo I, Civetta L, Dallai L, Moretti R, Orsi G, Andria M, Trecalli A (2013) Mantle and crustal processes in the magmatism of the Campania region: inferences from mineralogy, geochemistry, and Sr–Nd–O isotopes of young hybrid volcanics of the Ischia island (South Italy). Contrib Mineral Petr 165:1173–1194. https://doi.org/10.1007/s00410-013-0853-x

Davidson JP, McMillan NJ, Moorbath S, Worner G (1990) The Nevados de Payachata volcanic region (18S, 69 W, N. Chile) II. Evidence for widespread crustal involvement in Andean magmatism. Contrib Mineral Petr 105:412–432

Davidson JP, Tepley FJ III (1997) Recharge in volcanic systems; evidence from isotopic profiles of phe-nocrysts. Science 275:826–829

Davidson JP, Tepley FJ III, Knesel KM (1998) Crystal isotopic stratigraphy; a method for constraining magma differentiation pathways. Eos 79:185–189

De Campos CP, Dingwell DB, Fehr KT (2004) Decoupled convection cells from mixing experiments with alkaline melts from Phlegrean Fields. Chem Geol 213:227–251

De Campos CP, Dingwell DB, Perugini D, Civetta L, Fehr TK (2008) Heterogeneities in magma chambers: insight from the behaviour of major and minor elements during mixing experiments with natural alkaline melts. Chem Geol 256:131–145

De Campos CP, Perugini D, Ertel-Ingrisch W, Dingwell DB, Poli G (2011) Enhancement of Magma mixing efficiency by Chaotic Dynamics: an experimental study. Contrib Mineral Petr 161:863–881

Deino AL, Orsi G, de Vita S, Piochi M (2004) The age of the Neapolitan Yellow Tuff caldera-forming eruption (Campi Flegrei caldera—Italy) assessed by $^{40}$Ar/$^{39}$Ar dating method. J Volcanol Geotherm Res 133: 157–170

De Rosa R, Donato P, Ventura G (2002) Fractal analysis of min- gled/mixed magmas: an example from the Upper Pollara eruption (Salina Island, Southern Tyrrhenian Sea, Italy). Lithos 65:299–311

De Rosa R, Mazzuoli R, Ventura G (1996) Relationships between deformation and mixing processes in lava flows: a case study from Salina (Aeolian Islands, Tyrrhenian Sea). Bull Volcanol 58:286–297

de Vita S, Orsi G, Civetta L, Carandente A, D'Antonio M, Deino A, di Cesare T, Di Vito MA, Fisher RV, Isaia R, Marotta E, Necco A, Ort M, Pappalardo L, Piochi M, Southon J (1999) The Agnano-Monte Spina eruption (4100 years BP) in the restless Campi Flegrei caldera (Italy). J Volcanol Geotherm Res 91:269–301

Del Gaudio C, Aquino I, Ricciardi GP, Ricco C, Scan-done R (2010) Unrest episodes at Campi Flegrei: a reconstruction of vertical ground movements during 1905–2009. J Volcanol Geotherm Res 195:48–56

Di Muro A, Metrich N, Vergani D, Rosi M, Armienti P, Fourgeroux T, Deloule E, Arienzo I, Civetta L (2014) The shallow plumbing system of Piton de la Fournaise volcano (La Reunion island, Indian Ocean) revealed by the major 2007 caldera forming eruption. J Petrol
55(7):1287-1315. https://doi.org/10.1093/petrology/egu025

Di Renzo V, Arienzo I, Civetta L, D'Antonio M, Tonerini S, Di Vito MA, Orsi G (2011) The magmatic feeding system of the Campi Flegrei caldera: architecture and temporal evolution. Chem Geol 281: 227-241

Di Vito MA, Arienzo I, Braia G, Civetta L, D'Antonio M, Orsi G (2011) The Averno 2 fissure eruption: a recent small size explosive event at the Campi Flegrei caldera. Bull Volcanol 73:295-320. https://doi.org/10.1007/s00445-010-0417-0

Didier J, Barbarin B (1991) Enclaves and Granite Petrology, vol 13. Elsevier, Amsterdam

Eichelberger JC (1975) Origin of andesite and dacite: evidence of mixing at Glass Mountain in California and other circum-Pacific volcanoes. Geol Soc Amer Bull 86:1381-1391

Eichelberger JC (1978) Andesitic volcanism and crustal evolution. Nature 275:21-27

Eichelberger JC (1980) Vesiculation of mafic magma during replenishment of silicic magma reservoirs. Nature 288:446-450

Fedele L, Scarpati C, Lanphere M, Melluso L, Morra V, Perrotta A, Ricci G (2008) The Breccia Museo formation, Campi Flegrei, southern Italy: geochronology, chemostratigraphy and relationship with the Campanian Ignimbrite eruption. Bull Volcanol 70:1189-1219

Fedele L, Zanetti A, Morra V, Lustrino M, Melluso L, Vannucci R (2009) Clinopyroxene/liquid trace element partitioning in natural trachyte-trachyphonolite systems: insights from Campi Flegrei (southern Italy). Contrib Mineral Petr 158:337-356. https://doi.org/10.1007/s00410-009-0386-5

Flinders J, Clemens JD (1996) Non-linear dynamics, chaos, complexity and enclaves in granitoid magmas. Trans R Soc Edinburgh: Earth Sci 87:225-232

Font L, Davidson JP, Pearson DG, Nowell GM, Jer-ram DA, Ottley CJ (2008) Sr and Pb isotopic micro-analysis of plagioclase crystals from Skye lavas: an insight into open-system processes in a flood basalt province. J Petrol 49(8):1449-1471

Fourcade S, Allegre CJ (1981) Trace element behaviour in granite genesis: a case study the calc-alkaline plutonic association from the Querigut Complex (Pyrenees France). Contrib Mineral Petr 76:177-195

Francalanci L, Avanzinelli R, Nardini I, Tiepolo M, Davidson JP, Vannucci R (2012) Crystal recycling in the steady-state system of the active Stromboli volcano: a 2.5-ka story inferred from in situ Sr-isotopic and trace element data. Contrib Mineral Petr 163:109-131

Grasset O, Albaredé F (1994) Hybridisation of mingling magmas with different densities. Earth Planet Sci Lett 121:327-332

Hibbard MJ (1981) The magma mixing origin of mantled feldspar. Contrib Mineral Petrol 76:158-170

Hibbard MJ (1995) Petrography to petrogenesis. MacMil-lan Publishing Company, London

Jolis EM, Freda C, Troll VR, Deegan FM, Blythe LS, McLeod CL, Davidson JP (2013) Experimental simulation of magma-carbonate interaction beneath Mt. Vesuvius, Italy. Contrib Mineral Petr 166:1335-1353

Kent AJR, Darr C, Koleszar AM, Salisbury MJ, Cooper KM (2010) Preferential eruption of andesitic magmas through recharge filtering. Nat Geosci 3: 631-636

Kinman WS, Neal CR, Davidson JP, Font L (2009) The dynamics of Kerguelen Plateau magma evolution: new insights from major element, trace element and Sr isotopic microanalysis of plagioclase hosted in Elan Bank basalts. Chem Geol 264:247-265

Knesel KM, Davidson JP, Duffield WA (1999) Evolution of silicic magma through assimilation and subsequent recharge: evidence from Sr isotopics in sanidine phenocrysts Taylor Creek Rhyolites. J Petrol 40:773-786

Kress V (1997) Magma mixing as a source for Pinatubo sulphur. Nature 389:591-593

Lesher CE (1990) Decoupling of chemical and isotopic exchange during magma mixing. Nature 344:235-237

Longo A, Papale P, Vassalli M, Saccorotti G, Montagna CP, Cassioli A, Giudice S, Boschi E (2012) Magma convection and mixing dynamics as a source of Ultra-Long-Period oscillations. Bull Volcanol 74 (4):873-880. doi: https://doi.org/10.1007/s00445-011-0570-0, https://doi.org/10.1007/s00445-011

Martin VM, Morgan DJ, Jerram DA, Caddick MJ, Prior DJ, Davidson JP (2008) Bang! Month-scale eruption triggering at Santorini Volcano. Science 321:1178

Melluso L, De' Gennaro R, Fedele L, Franciosi L, Morra V (2012) Evidence of crystallization in residual, Cl-F-rich, agpaitic, trachyphonolitic magmas and primitive Mg-rich basalt-trachyphonolite interaction in the lava domes of the Phlegrean Fields (Italy). Geol Mag 149(3):532-550

Montagna CP, Papale P, Longo A (2015) Timescales of mingling in shallow magmatic reservoir. In: Caricchi L, Blundy JD (eds) Chemical, Physical and Temporal Evolution of Volcanic Systems. Geological Society, London, Special Publications 422, pp 131-140

Morgavi D, Perugini D, De Campos CP, Ertl-Ingrisch W, Dingwell DB (2013a) Time evolution of chemical exchanges during mixing of rhyolitic and basaltic melts. Contrib Mineral Petr 166(2):615-638

Morgavi D, Perugini D, De Campos CP, Ertl-Ingrisch W, Dingwell DB (2013b) Morphochemistry of patterns produced by mixing of rhyolitic and basaltic melts. J Volcanol Geotherm Res 253:87-96

Morgavi D, Perugini D, De Campos CP, Ertl-Ingrisch W, Lavallee Y, Morgan L, Dingwell DB (2013c) Interactions Between Rhyolitic and Basaltic Melts Unraveled by Chaotic Magma Mixing Experiments. Chem Geol 346:119-212

Morgavi D, Petrelli M, Vetere FP, Gonzalez-Gartia D, Perugini D (2015) High temperature apparatus for chaotic mixing of natural silicate melts. Rev Sci Instrum 86(10):105108
Morgavi D, Arzilli F, Pritchard C, Perugini D, Mancini L, Larson P, Dingwell BD (2016) The Grizzly Lake complex (Yellowstone Volcano, USA): mixing between basalt and rhyolite unraveled by microanalysis and X-ray microtomography. Lithos 260:457–474

Murphy MD (1998) The role of Magma Mixing in triggering the current eruption at the Soufriere Hills Volcano, Montserrat, West Indies. Geophys Res Lett 25:3433–3436

Laumonier M, Scaillet B, Pichavant M, Champallier R, Andujar J, Arbaret L (2014) On the conditions of magma mixing and its bearing on andesite production in the crust. Nature communications, volume 5

Papale P, Moretti R, Barbato D (2006) The compositional dependence of the saturation surface of H₂O + CO₂ fluids in silicate melts. Chem Geol 229(1–3):78–95

Pappalardo L, Piochi M, D'Antonio M, Civetta L, Petrini R (2002) Evidence for multistage magmatic evolution during the past 60 kyr at Campi Flegrei (Italy) deduced from Sr, Nd and Pb isotopic data. J Petrol 43(8):1415–1434

Perugini D, De Campos CP, Dingwell DB, Petrelli M, Poli G (2008) Trace element mobility during magma mixing: preliminary experimental results. Chem Geol 256:146–157

Perugini D, de Campos CP, Petrelli M, Dingwell DB (2015) Concentration variance decay during magma mixing: a volcanic chronometer. Sci Rep 5:14225

Perugini D, Petrelli M, Poli G (2006) Diffusive fractionation of trace elements by chaotic mixing of magmas. Earth Planet Sci Lett 243:669–680

Perugini D, Petrelli M, Poli G, De Campos C, Dingwell DB (2010) Time-scales of recent Phlegrean fields eruptions inferred from the application of a ‘Diffusive Fractionation’ model of trace elements. Bull Volcanol 72:431–447

Perugini D, Poli G (2005) Viscous fingering during replenishment of felsic magma chambers by continuous inputs of mafic magmas: field evidence and fluid-mechanics experiments. Geology 33:5–8

Perugini D, Poli G (2012) The mixing of magmas in plutonic and volcanic environments: analogies and differences. Lithos. https://doi.org/10.1016/j.lithos.2012.02.002

Perugini D, Poli G, Gatta G (2002) Analysis and simulation of Magma mixing processes in 3D. Lithos 65:313–330

Perugini D, Poli G, Mazzuoli R (2003) Chaotic advection, fractals and diffusion during mixing of Magmas: evidence from Lava Flows. J Volcanol Geotherm Res 124:255–279

Perugini D, Valentini L, Poli G (2007) Insights into Magma chamber processes from the analysis of size distribution of Enclaves in Lava Flows: a case study from Volcano Island (Southern Italy). J Volcanol Geotherm Res 166:193–203

Petrelli M, Perugini D, Poli G (2011) Transition to Chaos and implications for Time-scales of Magma Hybridization during mixing processes in Magma chambers. Lithos 125:211–220

Plail M, Barclay J, Humphreys MCS, Edmonds M, Herd RA, Christopher TE (2014) Characterization of mafic enclaves in the erupted products of Soufriere Hills Volcano, Montserrat, 2009–2010. In Wadge G, Robertson REA, Voight B, (eds) The eruption of Soufriere Hills Volcano, Montserrat from 2000 to 2010. Geological Society, London, Memoirs, 39, 343–360

Poli G, Tommasini S, Halliday AN (1996) Trace elements and isotopic exchange during acid-basic magma interaction processes. Trans Royal Soc Edinburgh: Earth Sci 87:225–232

Pritchard CJ, Larson PB, Spell TL, Tarbert KD (2013a) Eruption-triggered mixing of extra-caldera basalt and rhyolite complexes along the East Gallatin-Washburn fault zone, Yellowstone National Park, WY, USA. Lithos 175–176:163–177. doi:https://doi.org/10.1016/j.lithos.2013.04.022

Self S (1992) Krakatau revisited: the course of events and interpretation of the 1883 eruption. Geoj Lib 28:109–121

Sigmundsson F, Hreinsdottir S, Hooper A, Arnadottir T, Pedersen R, Roberts JM, Oskarsson N, Auriac A, Decriem J, Einarsson P, Geirsson H, Hensch M, Ofeigsson BG, Sturkell E, Sveinbjornsson H, Feigl LK (2010) Intrusion triggering of the 2010 Eyjafjallajokull explosive eruption. Nature 468:426–430

Slaby E, Gotze J, Worner G, Simon K, Wrzalik R, Smigielski M (2010) K-feldspar phenocrysts in micro-granular magmatic enclaves: a cathodoluminescence and geochemical study of crystal growth as a marker of magma mingling dynamics. Lithos 105:85–97

Smith JV (2000) Structures on interfaces of mingled magmas, Stewart Island, New Zealand. J Struct Geol 22:123–133

Sparks SRJ, Marshall LA (1986) Thermal and mechanical constraints on mixing between mafic and silicic magmas. J Volcanol Geotherm Res 29:99–124

Sparks SRJ, Sigurdsson H, Wilson L (1977) Magma mixing: a mechanism for triggering acid explosive eruptions. Nature 267:315–318

Tonarini S, D'Antonio M, Di Vito MA, Orsi G, Carandente A (2009) Geochemical and isotopical (B, Sr, Nd) evidence for mixing and mingling processes in the magmatic system feeding the Astroni volcano (4.1–3.8 ka) within the Campi Flegrei caldera (South Italy). Lithos 107:135–151

Tonarini S, Leeman WP, Civetta L, D'Antonio M, Ferrara G, Necco A (2004) B/Nb and 5¹¹B systematics in the Phlegrean Volcanic District (PVD). J Volcanol Geotherm Res 113:123–139

Ventura G (1998) Kinematic significance of mingling-rolling structures in lava flows: a case study from Porri volcano (Salina Southern Tyrrhenian Sea). Bull Volcanol 59:394–403

Wada K (1995) Fractal structure of heterogeneous ejecta from the Meakan volcano, eastern Hokkaido, Japan: implications for mixing mechanism in a volcanic conduit. J Volcanol Geotherm Res 66:69–79
Walker GPL, Skelhorn RR (1966) Some associations of
acid and basic igneous rocks. Earth Sci Rev 2:9–109

Wallace G, Bergantz G (2002) Wavelet-based correlation (WBC) of crystal populations and magma mixing. *Earth Planet Sci Lett* 202:133–145

Wiebe RA (1994) Silicic magma chambers as traps for basaltic magmas: the Cadillac mountain intrusive complex, Mount Desert island, Maine. J Geol 102:423–427

Wilcox RE (1944) Rhyolite-basalt complex of Gardiner River, Yellowstone Park, Wyoming. Geol Soc Am Bull 55:047–1080

Wilcox RE (1999) The idea of Magma Mixing: history of a struggle for Acceptance. J Geol 107(4):421–432

Williams-Jones G, Rymer H (2002) Detecting volcanic eruption precursos: a new method using gravity and deformation measurements. J Volcanol Geotherm Res 113:379–389

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.
<figure><img src="image_50.png" /></figure>

# Gases as Precursory Signals: Experimental Simulations, New Concepts and Models of Magma Degassing

M. Pichavant, N. Le Gall and B. Scaillet

## Abstract

Volatile release during magma ascent in volcanic conduits (magma degassing) forms the basis for using volcanic gases as precursory signals. Recent high temperature high pressure experimental simulations have yielded results that challenge key assumptions related to magma degassing and are important for the interpretation of glass inclusion and gas data and for using volcanic gas as precursory signals. The experimental data show that, for ascent rates expected in natural systems, pure H₂O basaltic melts will evolve mostly close to equilibrium when decompressed from 200 to 25 MPa. In the same way, degassing of H₂O-S species evolves at near equilibrium, although this conclusion is limited by the number of S solubility data available for basaltic melts. However, degassing of CO₂ is anomalous in all studies, whether performed on basaltic or rhyolitic melts. CO₂ stays concentrated in the melt at levels far exceeding solubilities. The anomalous behaviour of CO₂, when associated with near equilibrium H₂O losses, yields post-decompression glasses with CO₂ concentrations systematically higher than equilibrium degassing curves. Therefore, there is strong

experimental support for disequilibrium degassing during ascent of CO₂-bearing magmas. The existence of volatile concentration gradients around nucleated gas bubbles suggests that degassing is controlled by the respective mobilities (diffusivities) of volatiles within the melt. The recently formulated diffusive fractionation model reproduces the main characteristics, especially the volatile concentrations, of experimental glasses. The model also shows that the gas phase is more H₂O-rich than expected at equilibrium because CO₂ transfer toward the gas phase is hampered by its retention within the melt. However, only integrated gas compositions are calculated. Similarly, only bulk experimental fluid compositions are determined in recent experiments. Thus, constraints on the local gas phase are becoming necessary for the application to volcanoes. This stresses the need for the direct analysis of gas bubbles nucleated in decompression experiments. Pre-eruptive changes in volcanic CO₂/SO₂ and H₂O/CO₂ gas ratios are interpreted to reflect different pressures of gas-melt segregation in the conduit, an approach that assume gas-melt equilibrium. However, if disequilibrium magma degassing is accepted, the use of volatile saturation codes is no longer possible and caution must be exercised with the application of local equilibrium to volcanic gases. Future developments in the interpretation of gas data require progress

M. Pichavant (✉) · N. Le Gall · B. Scaillet
CNRS, Orléans, France
e-mail: michel.pichavant@cnrs-orleans.fr
from both sides, experimental and volcanolog-
ical. One priority is to reduce the gap in scales
between experiments and gas measurements.

**Keywords**

Magma ascent · Decompression experiments
Disequilibrium degassing · Volatile diffusion
Gas compositions

**Extended Abstract**

Volcanic gases are one of the main tools to
monitor changes in the activity of volcanoes and
forecast their eruption. Magma ascent toward the
surface is associated with the exsolution of vola-
tiles initially dissolved in the melt, a process
designated as “magma degassing”. Classically,
the interpretation of volcanic gases relies on the
assumption that degassing takes place at equilib-
rium. However, several observations (CO₂ con-
tents of basaltic seafloor glasses, H₂O and CO₂
concentrations in glass inclusions, explosive
basaltic volcanism) do not fit easily in such a
model. Recently, decompression, ascent and
degassing of magmas in volcanic conduits have
been simulated by high temperature high pressure
experiments. Results from these simulations
stress the need to critically reconsider the whole
mechanism of degassing in basaltic but also rhy-
olitic magmas. The new experimental data show
that, at the decompression rates tested, pure H₂O
basaltic melts will evolve mostly close to equi-
librium when decompressed from 200 to 25 MPa.
In the same way, degassing of H₂O–S species
evolves at near equilibrium, although this con-
clusion is limited by the number of S solubility
data available for basaltic melts. Degassing of
CO₂ is anomalous in all studies, whether per-
formed on basaltic or rhyolitic melts. CO₂ stays
concentrated in the melt at levels far exceeding
solubilities. The anomalous behaviour of CO₂,
when associated with near equilibrium H₂O los-
ses, yields post-decompression glasses with CO₂
concentrations systematically higher than equi-
librium degassing curves. Therefore, there is
strong experimental support for disequilibrium

degassing during ascent of CO₂-bearing magmas. The existence of volatile concentration gradients around nucleated gas bubbles suggests that degassing is controlled by the respective mobilities (diffusivities) of volatiles within the melt. The contrasted diffusivities of dissolved volatile species (in particular H₂O and CO₂) selectively limit their transfer toward the gas phase for timescales typical for magma ascent. The diffusive fractionation model recently formulated reproduces the main characteristics, especially the volatile concentrations, of experimental glasses. It provides a framework to interpret the new experimental observations and the systematic deviations from equilibrium observed in CO₂-bearing systems, although coupling between volatile diffusion and vesiculation requires a more elaborate treatment. The model also shows that the gas phase is more H₂O-rich than expected at equilibrium because CO₂ transfer toward the gas phase is hampered by its retention within the melt. However, only integrated gas compositions are calculated. In the same way, only bulk experimental fluid compositions are determined in recent experiments. Since the gas phase is essential for the application to volcanoes, constraints on the local gas phase are becoming necessary. Compositions of gas bubbles in decompression experiments must be linked not only with pressure but also with volatile concentrations of local melts, specific degassing textures and mechanisms. As a way in this direction, local gas-melt equilibrium assumes that chemical equilibrium persists at the local scale, despite evidence for disequilibrium at larger scales. However, there are alternative ways to constrain the composition of nucleated gas bubbles, thus stressing the need for their direct analysis in decompression experiments. Pre-eruptive changes in CO₂/SO₂ and, in some cases, H₂O/CO₂ gas ratios observed at several basaltic volcanoes are generally interpreted to reflect different pressures of gas-melt segregation in the conduit, an approach that assume gas-melt equilibrium. However, if disequilibrium magma degassing is accepted, volatile saturation codes can no longer be directly used. Caution also must be exercised with the application of local gas-melt equilibrium to volcanic gases which are probably
closer to integrated rather than to local compositions. Future developments in the interpretation of gas data require progress from both sides, experimental and volcanological. One priority is to reduce the gap in scales between experiments and gas measurements to refine interpretations of gas compositions as unrest signals.

# 1 Magma Degassing and Volcanic Gases as Precursory Signals

Volcanic gases are one of the main tools used to monitor changes in the activity of volcanoes and forecast their eruption. This approach is rooted in the strong pressure dependence of the solubility of volatiles (mainly H₂O, CO₂, SO₂, H₂S, Cl) in silicate melts. Accordingly, magma ascent toward the surface is associated with the exsolution of volatiles initially dissolved in the melt, a process designated as “magma degassing”. The different volatiles have contrasted solubilities in silicate melts and, therefore, are expected to react differently to decompression. This forms the basis for using volcanic gas ratios to infer magma ascent and depth of gas segregation in volcanic conduits. For example, the sudden increase of gas CO₂/SO₂ ratio has been used as an indication for deep magma recharge at Stromboli (Aiuppa et al. 2010). At Soufriere Hills volcano (Montserrat), a correlation has been noted between gas HCl/SO₂ and the level of shallow activity as marked by the rate of lava extrusion and dome growth (Christopher et al. 2010; Edmonds et al. 2010).

Classically, the interpretation of volcanic gases relies on the assumption that degassing takes place at equilibrium. In the case of basaltic magmas, this assumption is supported by the high temperatures, low viscosities and high volatile diffusivities (Sparks et al. 1994). Vesiculation (i.e., the combined processes of bubble nucleation, growth and coalescence) is thought to be relatively easy in basaltic melts and degassing of basaltic magma is classically viewed as an equilibrium process. However, several observations do not fit easily in

such a model. They include (1) the existence of basaltic seafloor glasses often supersaturated in CO₂ (e.g., Aubaud et al. 2004), (2) the occurrence of glass inclusions with H₂O and CO₂ concentrations inconsistent with closed system equilibrium degassing (e.g., Metrich et al. 2010) and (3) the occurrence of explosive basaltic volcanism (e.g., Head and Wilson 2003) which implies sudden rather than gradual release of volatiles.

Recently, decompression and ascent of basaltic magmas in volcanic conduits has been simulated by high temperature high pressure petrological experiments. These simulations stress the need to critically reconsider the whole mechanism of degassing in basaltic but also more silicic magmas. In particular, the assumption of equilibrium degassing is now becoming increasingly challenged. This has major implications for the interpretation of glass inclusion and gas data and, more generally, for the use of volcanic gas as precursory signals. In this Chapter, first, the recent experimental simulations are reviewed. We show that they all demonstrate an anomalous behaviour for CO₂ which tends to stay dissolved within the melt at concentrations too high for equilibrium. Second, the diffusive fractionation model which has been proposed to account for the new experimental observations is described and critically discussed. Finally, the implications of disequilibrium degassing for experimental fluid compositions and the interpretation of volcanic gas data as precursory signals are explored.

# 2 Experimental Simulations

## 2.1 Basaltic Systems

Following early work on systems with only pure CO₂ (Lensky et al. 2006), decompression experiments on hydrous basaltic melts have been carried out recently by Pichavant et al. (2013) at 1150–1180 °C, for initial pressures of 200–250 MPa, final pressures of 100, 50 and 25 MPa and for decompression rates between ~1.5 down to 0.25 m/s. Melts from Stromboli, pre-synthesized to incorporate dissolved H₂O (2.7–3.8 wt%) and CO₂ (600–1300 ppm), were
used as starting materials. The experiments were of continuous decompression type, and both constant (one ramp) and variable (two ramps) decompression rates were imposed. Final melt H₂O concentrations were homogeneous and always close to equilibrium solubility values. In contrast, the rate of vesiculation was found to control the final melt CO₂ concentration. High vesicularity charges had glass CO₂ concentrations that follow theoretical equilibrium degassing paths whereas glasses from low vesicularity charges showed marked deviations from equilibrium, with CO₂ concentrations up to one order of magnitude higher than equilibrium solubilities (Fig. 1a). The experimental results were interpreted in light of the slower diffusivity of CO₂ relative to H₂O in basaltic melts.

Yoshimura (2015) decompressed a natural evolved basaltic melt containing dissolved H₂O and CO₂ at 1200 °C and between 1000 and 500 MPa. The short decompression duration of 10 min over this pressure interval simulates a very fast ascent rate (~32 m/s for a rock density of 2650 kg/m³). A vesiculated glass was produced and Fourier Transform Infrared Spectroscopy (FTIR) profiles revealed large CO₂ concentration gradients in the melt adjacent to gas bubbles. In contrast, the melt H₂O content was almost constant throughout the sample. The glass volatile concentration data cover a near vertical trend in the H₂O–CO₂ diagram (Fig. 1a).

Le Gall and Pichavant (2016a) extended the decompression experiments performed by Pichavant et al. (2013), using essentially the same procedures and materials. Three starting volatile compositions were investigated: series #1 (4.91 wt% H₂O, no CO₂), series #2 (2.41 ± 0.04 wt% H₂O, 973 ± 63 ppm CO₂) and series #3 (0.98 ± 0.16 wt% H₂O, 872 ± 45 ppm CO₂). The volatile-bearing glasses were synthesized at 1200 °C and 200 MPa, then continuously decompressed at a fast decompression rate of 3 m/s in the pressure range 150–25 MPa and then rapidly quenched. Post-decompression glasses were characterized

texturally by X-ray microtomography. Volatile equilibrium was reached or approached during decompression in all series #1 melts with just water. In contrast, disequilibrium degassing occurred systematically in series #2 and #3 melts which retained elevated CO₂ concentrations (Fig. 1a). In similar experiments performed on the same three glass series but at a slower decompression rate of 1.5 m/s, Le Gall and Pichavant (2016b) found that series #1 (CO₂-free) melts followed equilibrium degassing until 100 MPa final pressure (Pfin). But at both 60 and 50 MPa Pfin, a slight H₂O-supersaturation was recognized, associated with a second bubble nucleation event that occurred at 25 MPa. In comparison, in series #2 and #3 (CO₂-bearing) melts, disequilibrium degassing was systematic, glasses retaining high non-equilibrium CO₂ concentrations (Fig. 1a).

The behavior of H₂O–, CO₂– and S-bearing basaltic melts during decompression was investigated by Le Gall et al. (2015a). Stromboli melts with 2.72 ± 0.02 wt% H₂O, 1291 ± 85 ppm CO₂ and 1535 ± 369 ppm S were synthesized at 1200 °C and 200 MPa and then decompressed to final pressures (Pfin) ranging from 150 to 25 MPa, followed by rapid quenching. The continuous decompressions were conducted at rates of 1.5 and 3 m/s. During decompression, S (and H₂O) were lost slightly more from the melt than expected from equilibrium degassing models, whilst significant CO₂ was retained at elevated concentrations in the melt (Fig. 1a). It was found that the degassing trend recorded by Stromboli glass inclusions could be closely reproduced by the experiments (Fig. 1b; Le Gall et al. 2015a). For andesitic melts, Fiege et al. (2014) observed that the fluid/melt partition coefficient for sulfur increases with the decompression rate. However, the influence of decompression rate on S degassing was marked only for oxidizing conditions, corresponding to sulfate as the only S species, thus making necessary to consider the different behaviour of S²⁻ and S⁶⁺ during degassing.
<figure><img src="image_51.png" /><figcaption>**Fig. 1** H₂O–CO₂ glass concentration diagrams for a basaltic decompression experiments, b Stromboli glass inclusions and c rhyolitic decompression experiments. In a, b and c, light grey curves are isobars labelled with pressure in MPa. In a, fields for post-decompression glasses are distinguished with (1) referring to Pichavant et al. (2013), (2) to Yoshimura (2015), (3) to Le Gall and Pichavant (2016a, b) and (4) to Le Gall et al. (2015a). Black curves are closed-system equilibrium degassing trajectories redrawn from the original figures (Pichavant</figcaption></figure>

et al. 2013; Yoshimura 2015; Le Gall and Pichavant 2016a, b; Le Gall et al. 2015a). In b, the glass inclusion data are from Metrich et al. (2010). MI: glass inclusions, Emb: embayments. In c, the glass data field and the bold theoretical equilibrium closed-system degassing curve are redrawn from Yoshimura (2015). Cₗₘᵢt and Cₗₙₑw are the composition of the pre-decompression melt and of the estimated post-decompression melt at the gas-melt interface, respectively (Yoshimura 2015; see also Figs. 5 and 6)
## 2.2 Rhyolitic Systems

Yoshimura (2015) decompressed a natural rhyolitic melt containing dissolved H₂O and CO₂ from 100 to 50 MPa at 800 °C. The duration of the decompression was 5000 s corresponding to a decompression rate of 0.38 m/s (for a density of 2650 kg/m³). FTIR analysis of the vesiculated glass sample showed CO₂ concentration gradients in the melt away from gas bubbles. In contrast, H₂O was found to be distributed homogeneously within the sample although H₂O concentrations decreased significantly relative to the pre-decompression melt showing it had reequilibrated. On the H₂O–CO₂ diagram, the glass volatile concentrations define a near vertical array located left of the theoretical equilibrium degassing curve (Fig. 1c).

## 2.3 Summary of Experimental Evidence

In all experiments above, melt vesiculation is the result of decompression, in most cases single-step (constant decompression rate) and, more rarely, multi-step (variable decompression rates, Pichavant et al. 2013). Vesiculation leads to the generation of a gas (or fluid) phase. Volatiles partition between melt and gas, and volatile concentrations in post-decompression glasses evolve from those initially dissolved in pre-decompression glasses. The evaluation of equilibrium vs. disequilibrium degassing is performed by comparing volatile concentrations of post-decompression glasses with theoretical closed-system equilibrium degassing trajectories (Fig. 1). For pure H₂O melts, this equilibrium trajectory is calculated using the experimental solubility data of Lesne et al. (2011). For CO₂- and S-bearing systems, gas-melt equilibrium thermodynamic models (Newman and Lowenstern 2002; Papale et al. 2006; Burgisser et al. 2015) are used. Results show that pure H₂O basaltic systems evolve close to equilibrium when decompressed from 200 to 25 MPa with ascent rates of 1.5 and 3 m/s, although small levels of H₂O supersaturation are observed below 100 MPa (Le Gall and Pichavant 2016a, b). In the same way,

degassing of S species evolves at near equilibrium (Le Gall et al. 2015a) although the reference equilibrium model (Burgisser et al. 2015) is somewhat uncertain due to the limited number of S solubility data for calibration. Degassing of CO₂ is anomalous in all studies, whether performed on basaltic or rhyolitic melts (Pichavant et al. 2013; Yoshimura 2015; Le Gall and Pichavant 2016a, b; Le Gall et al. 2015a). CO₂ stays concentrated in the melt at concentrations far exceeding solubilities (Fig. 1). Except in the very fast basalt decompression experiment of Yoshimura (2015), the anomalous behaviour of CO₂ is associated with significant H₂O losses which results in post-decompression glasses plotting systematically left of theoretical equilibrium degassing trajectories in Fig. 1. We conclude that recent experimental studies strongly support the possibility of disequilibrium degassing, i.e., that ascending melts can keep volatile concentrations (particularly CO₂) significantly different from those expected from equilibrium modelling. The questions thus arise of (1) the mechanisms responsible for this disequilibrium behaviour and (2) of the consequences of disequilibrium degassing for the composition of the gas phase.

# 3 Modelling Disequilibrium Degassing

## 3.1 The Diffusive Fractionation Model

Both Pichavant et al. (2013) and Yoshimura (2015) observed decoupling between the behaviour of H₂O and CO₂ during experimental decompression and degassing. In both studies, CO₂ concentration gradients were found in post-decompression glasses, either around gas bubbles or near the gas-melt interface. In contrast, no such diffusion profiles were identified for H₂O, despite concentrations being lower (in most cases) than in pre-decompression glasses (Pichavant et al. 2013; Yoshimura 2015; Le Gall et al. 2015a; Le Gall and Pichavant 2016a, b).

Pichavant et al. (2013) suggested that two characteristic distances, the gas interface distance
(either the distance between two bubbles in the melt or the distance to the gas-melt interface) and the volatile diffusion distance (a function of respective diffusivities of volatiles in the melt) control the degassing process. Yoshimura (2015) quantitatively formulated a diffusive fractionation model to describe the ascent and degassing of volatile-bearing magmas. The reader is referred to this work for details about the calculations. The model is based on a diffusivity of CO₂ being one log unit lower than for H₂O (e.g., Zhang and Ni 2010). Decompression trajectories computed from the model are shown on Fig. 2 for different ascent rates, from 0.1, 1, 10, 100 to ∞ m/s. Although very high ascent rates (e.g., Peslier et al. 2015) are necessary for degassing trajectories to shift significantly left to the equilibrium reference curve, the modelling results qualitatively reproduce the main characteristics of experimental post-decompression glasses, i.e., the elevated CO₂ glass concentrations, the significant H₂O losses and the melt concentration trends in H₂O–CO₂ diagrams (Fig. 1).

Yoshimura (2015) emphasized the relative simplicity of his model. For example, bubble growth was not considered as in other more elaborated theoretical treatments (e.g., Gonnermann and Manga 2005). Rather than continuously varying boundary (gas-melt) interface volatile concentrations and bubble-bubble distances as in a natural ascending magma, the calculations were performed step-by-step (i.e., at different pressures) along the decompression ramp, with fixed boundary concentrations and bubble-bubble distance (Yoshimura 2015). It is also important to note that the volatile concentrations on Fig. 2 correspond to averages computed by integrating the concentrations in the melt along the diffusion profiles (*distance integrated compositions*).

Gas phase compositions were calculated by mass balance using the initial volatile concentrations and the average volatile concentrations left in the melt after decompression and degassing (Yoshimura 2015). Results are shown on Fig. 3 and they correspond to compositions

<figure><img src="image_52.png" /><figcaption>**Fig. 2** H₂O and CO₂ melt volatile concentrations computed with the diffusive fractionation model for different decompression/ascent rates (v from 0.1 to ∞, in m/s). Isobars (light curves) are labelled with pressure in MPa.</figcaption></figure>

The heavy curve labelled “eq” is the equilibrium closed-system degassing trajectory as calculated by Yoshimura (2015). Figure redrawn from Yoshimura (2015). See text for details
<figure><img src="image_53.png" /><figcaption>**Fig. 3** Variations of the gas H₂O/CO₂ molar ratio (integrated compositions) with pressure computed from the diffusive fractionation model using different decompression/ascent rates (v from 0.1 to 1000, in m/s).</figcaption></figure>

The heavy curve labelled “eq” is the equilibrium closed-system degassing trajectory. Figure redrawn from Yoshimura 2015. See text for details

integrated along decompression (*pressure integrated compositions*). These compositions are more H₂O-rich (higher H₂O/CO₂ ratios) than gases generated under equilibrium degassing.

## 3.2 Coupling Between Diffusion and Vesiculation

Coupling between volatile diffusion and vesiculation is a necessity in diffusive degassing models because vesiculation defines the density of bubbles, their sizes and the distances between them (e.g., Pichavant et al. 2013; Le Gall et al. 2016a, b). This issue was addressed by Yoshimura (2015), although in a relatively simplified manner. The distance between bubbles was defined as being a function of only two variables, the distance between bubbles at the bottom of the decompression column (arbitrary value) and the vesicularity. Vesicularity must change along with decompression and degassing. So, the

vesicularity term should embody the textural variations associated with magma ascent. In the model of Yoshimura (2015), the vesicularity was computed from the amount of volatiles exsolved upon decompression, using an equation of state for H₂O–CO₂ gas mixtures to calculate the density of the gas and assuming a constant density for the melt. In so doing, it is apparent that only a vesicularity corresponding to equilibrium degassing is considered. Thus, for a given initial bubble-bubble distance, the distance between bubbles in the decompression column depends only on the equilibrium vesicularity. Degassing trajectories (Fig. 2) and integrated gas compositions (Fig. 3) were calculated on this basis.

For comparison, experimental vesicularities, bubble diameters and bubble number densities are shown on Fig. 4 for three series of basaltic melts decompressed from 200 to 25 MPa final pressure (Pfin) at 3 m/s (Le Gall and Pichavant 2016a). Systematic variations within the three glass series are observed depending on Pfin. In
most cases, the vesicularity data plot intermediate between the two equilibrium vesicularity curves, which were computed in a similar way than Yoshimura (2015) but only for two end-member cases corresponding to pure H₂O and pure CO₂ gas. The vesicularity data for the series #1 melts

(with pure H₂O) are in general much lower than the theoretical vesicularities calculated for pure H₂O gas. The data also show large changes in bubble sizes and bubble number densities that do not directly correlate with vesicularity. Le Gall and Pichavant (2016a) emphasized that

<caption>**Fig. 4** Textural data for post-decompression experimental glasses plotted as a function of final pressure (P<sub>fin</sub>) and comparison with data for natural basaltic pumices (Stromboli, Masaya). Vesicularities (**a**), bubble diameters (**b**) and bubble number densities (**c**) for three series of basaltic melts decompressed from 200 to 150, 100, 50 and to 25 MPa P<sub>fin</sub> at 3 m/s (Le Gall and Pichavant 2016a). Pre-decompression melt concentrations, series #1: 4.91 wt% H₂O (no CO₂), series #2: 2.41 ± 0.04 wt% H₂O, 973 ± 63 ppm CO₂ and series #3: 0.98 ± 0.16 wt% H₂O, 872 ± 45 ppm CO₂. Note that, for charge S + S38#1 which was partially fragmented, the vesicularity (**a**) and BND (**c**) data concern the unfragmented part. The bubble diameter data (**b**) are for both the unfragmented (black symbol) and the fragmented (minimum and maximum values) parts.</caption>

Figure redrawn from Le Gall and Pichavant (2016a)

<figure><img src="image_54.png" /></figure>
degassing textures result from several processes including bubble nucleation, growth, coalescence, plus buoyancy-driven bubble migration. We conclude that, although the diffusive fractionation model of Yoshimura (2015) provides a basis for coupling volatile diffusion calculations and vesiculation processes, more work is needed to incorporate the complex textural changes associated with ascent of volatile-bearing melts.

# 4 Implications for Gas Phase Compositions

## 4.1 Available Data and Models

Despite the limitations noted above, the diffusive fractionation model provides a framework to interpret the experimental observations and the systematic deviations from equilibrium degassing observed in CO₂-bearing systems. However, it should be emphasized that the model uses analytical data (glass volatile concentrations) and physicochemical properties (volatile diffusivities) related only to the melt phase. The question arises of the consequences of disequilibrium degassing for the gas phase composition. It is worth remembering here that the precursory signals come from gas data.

In the decompression experiments summarized above, the gas phase has not been chemically analysed although some mass balance calculations were performed to estimate the composition of the gas phase in the H₂O–, CO₂– and S-bearing experiments of Le Gall et al. (2015a). However, it is emphasized that, with this method, only bulk experimental gas compositions are provided (*charge and pressure integrated compositions*). No information is available on the composition of individual bubbles generated during decompression. The gas calculations performed by Yoshimura (2015) also use a similar mass balance approach, i.e., pressure integrated fluid compositions are given. However, the local gas at the gas-melt interface has an equilibrium composition (local gas-melt equilibrium). The differences between the disequilibrium (calculated with the model) and the equilibrium

(calculated assuming equilibrium degassing) gases (Fig. 3) is the consequence of CO₂ degassing being hampered by its retention within the melt. Therefore, disequilibrium is evidenced in the compositions of the pressure integrated fluids.

## 4.2 Composition of Gas Bubbles

The experiments and the diffusive fractionation model show that melt and gas both evolve under disequilibrium during magma ascent and degassing. For the melt, this conclusion is based either on volatile concentration measurements in glass at some distance of the gas/melt interface (Pichavant et al. 2013; Yoshimura 2015; Le Gall and Pichavant 2015a, 2016a, b) or on average concentrations calculated by integration along diffusion profiles (Yoshimura 2015). For the gas, constraints are available only on integrated compositions (Le Gall et al. 2015a; Yoshimura 2015). Since the gas phase is essential for the application to volcanoes, and given the interpretations proposed for the melt phase, constraints on the gas phase composition at smaller scales are becoming necessary. This requires linking compositions of gas bubbles in decompression experiments not only with pressure but also with volatile concentrations of local melts as well as with degassing textures and mechanisms.

As a way toward this direction, local gas-melt equilibrium can be assumed. This implies that chemical equilibrium persists locally between gas and melt, despite evidence for disequilibrium at larger scales. Therefore, the volatile compositions of melt and gas at the interface are defined by equilibrium partitioning of volatiles between these two phases (e.g., Dixon and Stolper 1995). To illustrate this concept, a schematic representation of the gas-melt interface for a H₂O– and CO₂-bearing melt decompressed isothermally from an initial (Pₘit) to a final (Pfin) pressure is shown on Fig. 5a. Initial volatile concentrations (Cₘit), together with the Pₘit and Pfin isobars and the equilibrium degassing trajectory are shown on the H₂O–CO₂ diagram of Fig. 5b. If local gas-melt equilibrium is assumed, the interface melt H₂O and CO₂ concentrations at Pfin (Cnew)
<figure><img src="image_55.png" /></figure>

**Fig. 5** Schematic illustration of local gas-melt equilibrium. a detail of the gas-melt interface region in a concentration ($C_i$) versus distance diagram where $C_i$ refers to the volatile concentration in the melt. The gas bubble is on the left. The two curves are melt volatile concentration profiles for H₂O (black) and CO₂ (grey) respectively, generated as a result of diffusion in the melt during decompression from $P_{init}$ to $P_{fin}$. $C_{init}$ give volatile concentrations of the pre-decompression melt, $C_{new}$ gas/melt interface volatile concentrations at $P_{fin}$ and $C_{sol}$ volatile solubilities at $P_{fin}$. Black lettering is used for

are equal to their solubilities ($C_{sol}$) at $P = P_{fin}$ (intersection of the equilibrium degassing curve with the $P_{fin}$ isobar, Fig. 5b). Note that diffusive fractionation generates H₂O and CO₂ concentration gradients within the melt (Fig. 5a), the range of possible melt compositions during decompression being represented by the dark grey domain in Fig. 5b. The interface melt is the only melt at equilibrium with the local gas at $P_{fin}$ which has a CO₂/H₂O corresponding to the fluid isopleth on Fig. 5b (e.g., Dixon and Stolper 1995). For a low pressure (e.g., 25 MPa), the local gas (e.g., a gas bubble nucleated at $P_{fin}$) is relatively H₂O-rich. In comparison, the pressure integrated gas assuming bulk *equilibrium* degassing from $P_{init}$ to $P_{fin}$ would be necessarily less H₂O-rich since most of the CO₂ must have been outgassed from the melt. This gas is less H₂O-rich than the pressure integrated gas produced by *disequilibrium* degassing from $P_{init}$ to $P_{fin}$ (Fig. 3). Thus, individual bubbles nucleated at $P_{fin}$ can have CO₂/H₂O different from the composition of integrated gases generated continuously during decompression.

H₂O and grey for CO₂. b H₂O–CO₂ diagram illustrating the evolution during decompression and degassing. The black bold curve is the equilibrium degassing trajectory. The two black lines are isobars labelled with initial ($P_{init}$) and final ($P_{fin}$) pressures along the decompression path. The dashed curve is the CO₂/H₂O isopleth passing through $C_{sol}$ and it defines the composition of the gas bubble in local equilibrium with the interface melt. The shaded domain gives the range of possible melt compositions generated upon decompression from $P_{init}$ down to $P_{fin}$

An alternative way to constrain the composition of gas bubbles is illustrated on Fig. 6. It starts from the observation that bubble nucleation is, from a kinetic point of view, an instantaneous process (e.g., Mourtada-Bonnefoi and Laporte 2002, 2004). Nucleation of a gas bubble draws volatiles from the local melt and the possibility that the initial CO₂/H₂O of the gas bubble is the same as the local melt should be considered. According to this hypothesis, represented schematically on Fig. 6, the local melt next to the nucleated bubble ($C_{new}$) is volatile-depleted but its CO₂/H₂O (r, Fig. 6a) is the same than the initial melt ($C_{init}$). Melt and gas bubble compositions are thus both located on a mixing line between $C_{init}$ and $C_{new}$ which passes through the origin of the H₂O–CO₂ diagram (Fig. 6b). The net result is the nucleation of individual gas bubbles more CO₂-rich than expected from local gas-melt equilibrium (Fig. 5b).

The previous discussion emphasizes the compositional variability of nucleated gas bubbles and the need for their direct analysis in decompression experiments. Comparison between the
<figure><img src="image_56.png" /></figure>

**Fig. 6** Hypothetical model for the composition of a gas bubble nucleated during decompression of H₂O–, CO₂-bearing basaltic melts. **a** detail of the gas-melt interface region in a concentration (Cᵢ) versus distance diagram where Cᵢ refers to the volatile concentration in the melt. The gas bubble is on the left. The two horizontal lines are melt volatile concentrations for H₂O (black) and CO₂ (grey) drawn as straight lines because bubble nucleation is an instantaneous event. The narrow zone near the interface is the melt region depleted in volatiles drawn to form the bubble. Cₙew are volatile concentrations in the depleted melt region. Notice that the H₂O/CO₂ ratio (r) is identical in both the depleted and non-depleted melt regions because during nucleation volatiles are drawn from the local melt and the initial CO₂/H₂O of the gas bubble is the same as the local melt. Cₘit give volatile concentrations of the

pre-decompression melt, C<sub>new</sub> gas/melt interface volatile concentrations left after bubble nucleation at P<sub>fin</sub>. Black lettering is used for H<sub>2</sub>O and grey for CO<sub>2</sub>. b H<sub>2</sub>O–CO<sub>2</sub> diagram illustrating the evolution during decompression and degassing. The black bold curve gives the schematic location of the theoretical equilibrium degassing trajectory. The two black lines are isobars labelled with initial (P<sub>init</sub>) and final (P<sub>fin</sub>) pressures along the decompression path. Melts produced as a result of decompression and bubble nucleation plot on the straight line joining C<sub>init</sub> and C<sub>new</sub>. This line passes through the origin of the diagram because both C<sub>init</sub> and C<sub>new</sub> have the same CO<sub>2</sub>/H<sub>2</sub>O ratio. Note that the location of C<sub>new</sub> along this line is arbitrary. The CO<sub>2</sub>/H<sub>2</sub>O ratio of the gas bubble (dashed line) is also the same as C<sub>init</sub> and C<sub>new</sub>. It is higher than the gas bubble controlled by local gas-melt equilibrium (Fig. 5).

CO₂/H₂O of nucleated bubbles and results of gas-melt volatile partitioning models (e.g., Dixon and Stolper 1995; Papale et al. 2006) would provide a crucial test of the local gas-melt equilibrium model (Fig. 5). If the nucleated bubbles prove to be CO₂-rich, then alternative models of control of gas composition would be supported (Fig. 6). One important aspect is that, in CO₂-bearing systems, bubble nucleation during decompression is continuous, occurring over a large pressure range (Le Gall and Pichavant 2016a, b). This is because CO₂-bearing melts are volatile-supersaturated (Fig. 1) and, so, the driving force for nucleation of new bubbles is always present. Thus, decompression of CO₂-bearing melts continuously leads to the nucleation of new bubbles, which increases the relevance of the hypothetical mechanism illustrated in Fig. 6.

# 5 Discussion and Perspectives for Gas Monitoring

## 5.1 Degassing Processes

Experimental simulations show that, for ascent rates expected in natural systems, equilibrium degassing occurs in pure H₂O melts. In contrast, results for CO₂-bearing melts conclusively demonstrate that degassing generates melt volatile concentrations out of equilibrium. The experimental database supporting this conclusion has recently expanded. It now includes basaltic and rhyolitic melts, and S-bearing as well as S-free systems. Several of those experimental decompression studies have been scaled to natural systems so that results are realistic and
applicable. The decompression experiments on Stromboli basalt cover ascent rates ranging from 0.25 to 3 m/s (Pichavant et al. 2013; Le Gall and Pichavant 2016a, b; Le Gall et al. 2015a), well in the range of current estimates for basaltic magmas (e.g., Rutherford 2008; Peslier et al. 2015). An ascent rate approximately 10 times faster was used by Yoshimura (2015). Therefore, in the case of basaltic magmas, equilibrium degassing should be viewed more as a reference situation rather than as a general mechanism. This is quite a change in paradigm which has major implications for how gas signals are interpreted.

One remaining issue concerns the role of crystals on bubble nucleation, heterogeneous rather than homogeneous. All decompression studies considered in this paper were performed on very crystal-poor, if not totally crystal-free, melts and bubble nucleation appears to be mostly homogeneous. Crystals present in experimental basaltic products include Fe-Ti oxides (Le Gall and Pichavant 2015b) and rare Fe sulphides (Le Gall et al. 2015a). Le Gall and Pichavant (2015b) have documented heterogeneous nucleation of bubbles on Fe-Ti oxide crystals (and also on Fe sulphides, Le Gall et al. 2015a). Recently, Shea (2017) has stressed the importance of magnetite as a key mineral phase promoting heterogeneous bubble nucleation in natural magmas. However, Fe oxide phenocrysts and sulphides are uncommonly present in amounts exceeding a few vol.% in natural magmas. This led Le Gall and Pichavant (2015b) to conclude that heterogeneous bubble nucleation is not an important mechanism in basaltic melts if driven by Fe oxides. Yet, heterogeneous nucleation on silicate phases is still an open question. For example, olivine, clinopyroxene and plagioclase are typical phenocrysts and microlites in Stromboli basalts (e.g., Pichavant et al. 2011). On the basis of limited textural evidence, Pichavant et al. (2013) ruled out the possibility of heterogeneous nucleation of gas bubbles on clinopyroxene and olivine crystals. However, additional investigations seem warranted to guarantee full applicability of the decompression experiments above.

Disequilibrium degassing, as documented in the experiments, is the consequence of the anomalous behaviour of CO₂. CO₂-super-saturated melts are systematically generated during decompression. The interpretation suggested by Pichavant et al. (2013) and quantitatively formulated by Yoshimura (2015) is that, because of its restricted diffusive mobility within the melt, CO₂ has limited access to the gas phase for timescales typical of magma ascent. However, our knowledge of volatile diffusivities in silicate melts is still very fragmentary. There are very few diffusivity data for H₂O and CO₂ on the same melt. S is another volatile which reputedly has a slow diffusivity in silicate melts. Yet, the behaviour of S during degassing differs from that of CO₂, although we are still short of S solubility data for basaltic melts (e.g., Lesne et al. 2015). Acquisition of fundamental data (especially volatile diffusivity and solubility data for basaltic melts) is needed for the elaboration of more detailed interpretations of the decompression experiments. Future magma ascent models should also incorporate the textural complexities associated with the vesiculation process.

## 5.2 Gases as Unrest Signals

Pre-eruptive changes in gas ratios have been observed at several basaltic volcanoes such as Stromboli (Burton et al. 2007; Aiuppa et al. 2010), Etna (Aiuppa et al. 2007) and Villarrica (Aiuppa et al. 2017) among others. Transition from passive degassing to more explosive paroxysmal eruption regimes is marked by temporal increases of the CO₂/SO₂ gas ratio in the volcanic plume. In some cases, the CO₂/SO₂ variations are correlated with a decrease of the H₂O/CO₂ gas ratio (e.g., Aiuppa et al. 2017). These variations in volcanic gas ratios have been generally interpreted to reflect different pressures of gas-melt segregation in the conduit, high CO₂/SO₂ (and low H₂O/CO₂) indicating deep conditions and low CO₂/SO₂ (and high H₂O/CO₂) shallow conditions (e.g., Edmonds 2008; Burton
et al. 2007; Allard 2010; Aiuppa et al. 2017). In this approach, the pressure-dependent evolution of the gas phase exsolved upon magma ascent and decompression is calculated by using volatile saturation codes (Newman and Lowenstern 2002; Moretti and Papale 2004; Papale et al. 2006; Burgisser et al. 2015). This implicitly assumes chemical equilibrium between gas and melt, an assumption which, as shown above, is now largely questioned. If disequilibrium magma degassing is accepted, then the consequences for the interpretation of gas signals need to be examined.

Firstly, one might argue that gas-melt equilibrium can persist at local scale, despite disequilibrium at larger scales. Thus, volatile saturation codes could still be used and applied to local gas and melt compositions, for example to model the composition of unconnected bubbles nucleated within the melt. In contrast, volcanic gases necessarily require, to be sampled, that the magma is permeable and, so, that the gas phase is connected. It is quite possible that the gases sampled are mixtures of different components, either integrated from several discrete degassing events along ascent or issued from different parts of the plumbing system. Therefore, volcanic gases are probably more representative of integrated compositions as discussed above than to compositions of local gases. We have shown previously that individual bubbles with compositions defined by local gas-melt equilibrium at a given pressure (Fig. 5) can have CO₂/H₂O different from integrated gases generated continuously during decompression (Fig. 3). We conclude to the limited applicability of local gas melt equilibrium to interpret volcanic gas ratios.

Secondly, disequilibrium gas-melt degassing due to CO₂ retention within the melt implies that CO₂/SO₂ and H₂O/CO₂ gas ratios can no longer be directly related to pressures of gas-melt segregation. Calculations using the diffusive fractionation model (Fig. 3) show that the pressure integrated gases have a higher H₂O/CO₂ (and also presumably a lower CO₂/SO₂ because CO₂ is retained

within the melt) than the same gases calculated assuming equilibrium with the melt (Fig. 3). This demonstrates the possibility of changing the gas ratios depending on the degassing mechanism (equilibrium vs. disequilibrium). It is worth emphasizing that disequilibrium degassing associated with CO₂ retention produces integrated fluids that are less, not more, CO₂-rich (Fig. 3).

The CO₂-rich gases observed on basaltic volcanoes have been generally attributed to deep-seated processes such as fluxing of CO₂ or arrival of CO₂-rich magmas (e.g., Aiuppa et al. 2010, 2017; Allard 2010). In contrast, the degassing mechanism of Fig. 6 (although it needs validation from direct analysis of gas bubbles in decompression experiments) allows CO₂-rich gas bubbles to be generated at low pressures. It also provides an example of how gas ratios can be changed at constant pressure depending on the degassing mechanism. The initially CO₂-rich bubbles (Fig. 6) will probably shift rapidly with time toward lower CO₂/H₂O because of preferential diffusion of H₂O from the melt. However, nucleation is a continuous process in CO₂-bearing basaltic melts (Le Gall and Pichavant 2016a, b) and reequilibration of previously nucleated bubbles by diffusion will be accompanied by the nucleation of new CO₂-rich bubbles.

We conclude that future developments in the interpretation of gas data require progress from both sides, experimental and volcanological. Some crucial experimental information at small scale is still missing such as the composition of individual gas bubbles nucleated in the decompression experiments and the influence of crystals on bubble nucleation. In parallel, at larger scales, the representativity and the significance of the gas phase sampled on active basaltic volcanoes needs to be better demonstrated, for example by combining gas measurements with detailed textural studies of eruption products. It is expected that future work will narrow the gap in scales between experiments and gas measurements to refine interpretations of gas compositions as unrest signals.
**Acknowledgements** This paper has benefited from discussions with P. Allard, C. Martel, N. Metrich, A. Bertagnini, R. Moretti, P. Papale and M. Pompilio, reviews by R. Brooker and F. Wadsworth and from editorial comments by B. Scheu. Discussion with S. Yoshimura was helpful. The VUELCO consortium provided a scientifically demanding and interdisciplinary forum for the elaboration of ideas developed in this study. The Ph.D. thesis of NLG was supported by the VUELCO project.

## References

Aiuppa A, Moretti R, Federico C, Giudice G, Gurrieri S, Liuzzo M, Papale P, Shinohara H, Valenza M (2007) Forecasting Etna eruptions by real-time observation of volcanic gas composition. Geology 35:1115–1118

Aiuppa A, Bertagnini A, Metrich N, Moretti R, Di Muro A (2010) A model of degassing for Stromboli volcano. Earth Planet Sci Lett 295:195–204

Aiuppa A, Bitetto M, Francofonte V, Velasquez G, Bucarey Parra C, Giudice G, Liuzzo M, Moretti R, Moussallam Y, Peters N, Tamburello G, Valderama OA, Curtis A (2017) A CO₂-gas precursor to the March 2015 Villarrica volcano eruption. Geochem Geophys Geosyst 18:2120–2132

Allard P (2010) A CO₂-rich gas trigger of explosive paroxysms at Stromboli basaltic volcano, Italy. J Volcanol Geotherm Res 189:363–374

Aubaud C, Pineau F, Jambon A, Javoy M (2004) Kinetic disequilibrium of C, He, Ar and carbon isotopes during degassing of mid-ocean ridge basalts. Earth Planet Sci Lett 222:391–406

Burgisser A, Alletti M, Scaillet B (2015) Simulating the behavior of volatiles belonging to the C-O-H-S system in silicate melts under magmatic conditions with the software D-Compress. Comput Geosci 79:1–14

Burton M, Allard P, La Spina A, Murè F (2007) Magmatic gas composition reveals the source depth of slug-driven strombolian explosive activity. Science 317:227–230

Christopher T, Edmonds M, Humphreys MCS, Herd RA (2010) Volcanic gas emissions from Soufrière Hills Volcano, Montserrat 1995–2009, with implications for mafic magma supply and degassing. Geophys Res Lett 37: L00E04. https://doi.org/10.1029/2009gl041325

Dixon JE, Stolper EM (1995) An experimental study of water and carbon dioxide solubilities in mid-ocean ridge basaltic liquids. Part II. Applications to degassing. J Petrol 36:1633–1646

Edmonds M (2008) New geochemical insights into volcanic degassing. Philos Trans R Soc A 366:4559–4579

Edmonds M, Aiuppa A, Humphreys M, Moretti R, Giudice G, Martin RS, Herd RA, Christopher T (2010) Excess volatiles supplied by mingling of mafic magma

at an andesite arc volcano. Geochem Geophys Geosyst 11:Q04005. https://doi.org/10.1029/2009GC002781

Fiege A, Behrens H, Holtz F, Adams F (2014) Kinetic vs. thermodynamic control of degassing of H₂O-S ± Cl-bearing andesitic melts. Geochim Cosmochim Acta 125:241–264

Gonnermann HM, Manga M (2005) Non-equilibrium magma degassing: results from modelling of the ca. 1340 AD eruption of Mono craters, California. Earth Planet Sci Lett 238:1–16

Head JW III, Wilson L (2003) Deep submarine pyroclastic eruptions: theory and predicted landforms and deposits. J Volc Geotherm Res 121:155–193

Le Gall N, Pichavant M (2015b) Heterogeneous bubble nucleation on Fe-Ti oxides in H₂O- and H₂O-CO₂-bearing basaltic melts. (In preparation)

Le Gall N, Pichavant M (2016a) Homogeneous bubble nucleation in H₂O- and H₂O-CO₂-bearing basaltic melts: results of high temperature decompression experiments. J Volcanol Geotherm Res 327:604–621

Le Gall N, Pichavant M (2016b) Effect of ascent rate on homogeneous bubble nucleation in the system basalt-H₂O-CO₂: Implications for Stromboli volcano. Am Mineral 101:1967–1985

Le Gall N, Pichavant M, Di Carlo I, Scaillet B (2015a) Sulfur partitioning between melt and fluid during degassing of ascending C-O-H-S-bearing basaltic magma: an experimental study. In preparation

Lensky NG, Niebo RW, Holloway JR, Lyakhovsky V, Navon O (2006) Bubble nucleation as a trigger for xenolith entrapment in mantle melts. Earth Planet Sci Lett 245:278–288

Lesne P, Scaillet B, Pichavant M, Iacono-Marziano G, Bény J-M (2011) The H₂O solubility of alkali basaltic melts: an experimental study. Contrib Mineral Petrol 162:133–151

Lesne P, Scaillet B, Pichavant M (2015) The solubility of sulphur in hydrous basaltic melts. Chem Geol 418:104–116

Metrich N, Bertagnini A, Di Muro A (2010) Conditions of magma storage, degassing and ascent at Stromboli: new insights into the volcanic plumbing system with inferences on the eruptive dynamics. J Petrol 51:603–626

Moretti R, Papale P (2004) On the oxidation state and volatile behaviour in multicomponent gas-melt equilibrium. Chem Geol 213:265–280

Mourtada-Bonnefoi CC, Laporte D (2002) Homogeneous bubble nucleation in rhyolitic magmas: an experimental study of the effect of H₂O and CO₂. J Geophys Res 107(B4):2066. https://doi.org/10.1029/2001jb000290

Mourtada-Bonnefoi CC, Laporte D (2004) Kinetics of bubble nucleation in a rhyolitic melt: an experimental study of the effect of ascent rate. Earth Planet Sci Lett 218:521–537

Newman S, Lowenstern JB (2002) VolatileCalc: a silicate melt-H₂O-CO₂ solution model written in Visual Basic for Excel. Comput Geosci 28:597–604
Papale P, Moretti R, Barbato D (2006) The compositional dependence of the saturation surface of H₂O + CO₂ fluids in silicate melts. *Chem Geol* 229:78–95

Peslier AH, Bizimis M, Matney M (2015) Water disequilibrium in olivines from Hawaiian peridotites: recent metasomatism, H diffusion and magma ascent rates. *Geochim Cosmochim Acta* 154:98–117

Pichavant M, Pompilio M, D'Oriano C, Di Carlo I (2011) The deep feeding system of Stromboli, Italy: insights from a primitive golden pumice. *Eur J Miner* 23:499–517

Pichavant M, Di Carlo I, Rotolo SG, Scaillet M, Burgiser A, Le Gall N, Martel C (2013) Generation of CO₂-rich melts during basalt magma ascent and degassing. *Contrib Mineral Petrol* 166:545–561

Rutherford MJ (2008) Magma ascent rates. In: Putirka KD, Tepley F (eds) Minerals, inclusions and volcanic processes. *Mineralogical Society of America Reviews in Mineralogy* vol 69, pp 241–271

Shea T (2017) Bubble nucleation in magmas: A dominantly heterogeneous process? *J Volcanol Geotherm Res* 343:155–170

Sparks RSJ, Barclay J, Jaupart C, Mader HM, Phillips JC (1994) Physical aspects of magmatic degassing I. Experimental and theoretical constraints on vesiculation. In: Carroll MR, Holloway JR (eds) Volatiles in Magmas. *Mineralogical Society of America Reviews in Mineralogy*, vol 30, pp 413–445

Yoshimura S (2015) Diffusive fractionation of H₂O and CO₂ during magma degassing. *Chem Geol* 411:172–181

Zhang Y, Ni H (2010) Diffusion of H, C, and O components in silicate melts. In: Zhang Y, Cherniak DJ (eds) Diffusion in minerals and melts. *Mineralogical Society of America Reviews in Mineralogy* vol 72, pp 171–225

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

<figure><img src="image_57.png" /></figure>
# Crystals, Bubbles and Melt: Critical Conduit Processes Revealed by Numerical Models

M. E. Thomas, J. W. Neuberg and A. S. D. Collinson

## Abstract

Understanding how magma moves within a conduit is an important question that is still poorly understood. In particular, estimation of the magma ascent rate is key for interpreting monitoring signals and therefore, predicting volcanic activity. This relies on understanding how strongly different magmatic processes occurring within the conduit control the ascent rate. These processes are controlled by changes in magmatic parameters such as the water content or temperature and understanding/ linking changes of such parameters to moni- toring data is an essential step in the use of these data as a predictive tool. The results presented here are from a suite of conduit flow models based on Soufrière Hills Volcano, Montserrat, that assesses the influence of individual model parameters. By systematically changing these parameters, the results indicate that changes in conduit diameter and excess pressure in the magma chamber are amongst the dominant controlling variables. However, the single most important parameter controlling variations in the magma ascent rate is the

volatile content. Therefore, understanding the processes controlling the volatile content within the conduit system and the outgassing of these volatiles is crucial to understanding and pre- dicting potential unrest or eruption scenarios.

## Keywords

Numerical modelling • Conduit processes
Low frequency earthquakes • Magma flow
Magma ascent rate

## 1 Introduction

A volcanic conduit provides the pathway for transport of magma and magmatic fluids within a volcano. It is possible to detect both this move- ment and the occurrence of conduit processes through geophysical monitoring techniques as discussed in Chap. 2 of this book. However, the extent to which changes in magma flow proper- ties affect the data recorded on volcanoes is not well understood. Is it possible that a small change in magma temperature or water content could alter the processes or flow within the conduit enough to be recorded by geophysical monitor- ing instruments or simple visual observation? What effect does the size of gas bubbles within

<footnote>M. E. Thomas (✉) · J. W. Neuberg
A. S. D. Collinson
School of Earth and Environment, Institute
of Geophysics and Tectonics, University of Leeds,
Leeds LS12 9JT, UK
e-mail: m.e.thomas@leeds.ac.uk</footnote>
the magma have on the overall flow dynamics, and how big do these changes need to be to alter the eruption style? These types of question are addressed within this chapter in an attempt to identify the crucial parameters that cause changes in observed volcanic behaviour.

We use conduit flow models to analyse the key input parameters that control magma flow properties, such as the magma water content, crystal content and conduit geometry, to assess their relative importance to the overall magma flow dynamics. A list of all input parameters is presented in Table 1 along with the range of values studied. We focus on evolved silicic magmatic systems because of the wealth of relevant monitoring information and previous numerical modelling attempts relating to

Soufrière Hills Volcano, Montserrat—a long lived andesitic dome forming eruption (Sparks et al. 2000; Wadge et al. 2014) and excellent natural laboratory. While these initial models are based on extrusive eruptions, the results of changing the model parameters have the potential to alter the eruption style to either more violent or gentle forms and it is noted that the underlying principles discussed here are applicable to other volcanic systems, including those exhibiting signs of unrest, yet to develop into a full blown eruption. If we can develop the findings presented in this chapter into the creation of threshold levels for recorded geophysical data, it may become possible to begin to predict when a volcano will evolve from a state of unrest to eruption.

<table><caption>**Table 1** Parameters used in the reference model and range of parameter variations</caption><thead><tr><th>Symbol/abbreviation</th><th>Variable</th><th>“Reference” model value</th><th>Range of modelled values</th></tr></thead><tbody><tr><td>–</td><td>The melt composition</td><td>Rhyolitic (&gt;71% SiO<sub>2</sub>) (Barclay et al. 1998)</td><td>See Table 2</td></tr><tr><td>b<sub>ni</sub></td><td>Bubble number density</td><td>10<sup>10</sup> m<sup>-3</sup> (Cluzel et al. 2008)</td><td>10<sup>7</sup>–10<sup>11</sup> m<sup>-3</sup></td></tr><tr><td>D<sub>TBL</sub></td><td>Thickness of thermal boundary layer over which T<sub>diff</sub> is lost</td><td>0.3 m (Collier and Neuberg 2006)</td><td>0.3–0.5 m</td></tr><tr><td>Γ</td><td>Bubble surface tension</td><td>0.06 N m<sup>-1</sup> (Lyakhovsky et al. 1996)</td><td>0.05–0.25 N m<sup>-1</sup></td></tr><tr><td>χ<sub>c</sub></td><td>Magma chamber crystal volume fraction</td><td>40% (Barclay et al. 1998)</td><td>40–50%</td></tr><tr><td>Ls</td><td>Slip length of brittle failure of melt</td><td>0.01</td><td>0.01–1.0 m</td></tr><tr><td>P<sub>e</sub></td><td>Excess chamber pressure above lithostatic</td><td>0 MPa</td><td>0–20 MPa</td></tr><tr><td>P<sub>top</sub></td><td>Pressure at conduit exit</td><td>0.09 MPa</td><td>0.09–4.5 MPa</td></tr><tr><td>ρ<sub>c</sub></td><td>Average density of crystal assemblage</td><td>2700 kg m<sup>-3</sup> (Burgisser et al. 2010)</td><td>2550–3200 kg m<sup>-3</sup></td></tr><tr><td>ρ<sub>m</sub></td><td>Density of pure melt</td><td>2380 kg m<sup>-3</sup> (Burgisser et al. 2010)</td><td>–</td></tr><tr><td>T</td><td>Magma temperature</td><td>1150 K (Devine et al. 2003)</td><td>1100–1150 K</td></tr><tr><td>T<sub>diff</sub></td><td>Amount of cooling at conduit wall</td><td>200 K (Collier and Neuberg, 2006)</td><td>100–200 K</td></tr><tr><td>τ<sub>s</sub></td><td>Melt shear strength</td><td>–</td><td>10<sup>5</sup>–10<sup>7</sup> Pa</td></tr><tr><td>W<sub>%</sub></td><td>Initial dissolved water content of magma</td><td>4.5 wt% (Barclay et al. 1998)</td><td>3–8 wt%</td></tr><tr><td>w, d, r</td><td>Variables that define the conduit shape and size</td><td>See Fig. 1</td><td>See Fig. 2</td></tr></tbody></table>
## 2 The Model

In order to assess the effect of altering the model parameters, a standard or “reference” model is defined. This reference model is based on data available in the literature that refers to Soufrière Hills Volcano, and is outlined in Fig. 1 and Table 1. The general dimensions of the modelled conduit, shown in Fig. 1 are inferred from geochemical and observational data from Soufrière Hills Volcano (Barclay et al. 1998; Sparks et al. 2000), placing minimum depth constraints of 5–6 km for the position of the magma chamber and width estimates of 30–50 m for the conduit.

### 2.1 Governing Equations

Conduit flow is computed with a finite element approach within the code COMSOL Multiphysics®, and modelled in an axial symmetric domain space through the compressible formulation of the Navier-Stokes equation:

$$ \rho \frac{\partial \mathbf{u}}{\partial t} + \rho \mathbf{u} \cdot \nabla \mathbf{u} = -\nabla p + \nabla \cdot \{\eta [\nabla \mathbf{u} + (\nabla \mathbf{u})^T] - \frac{2}{3} \eta [\nabla \cdot \mathbf{u}] \mathbf{I}\} + \mathbf{F} \quad (1) $$

and the continuity equation:

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0 \quad (2) $$

<figure><img src="image_58.png" /><figcaption>**Fig. 1** Cartoon of the modelled volcanic system. Bubble nucleation and brittle failure depth vary with the model parameters considered</figcaption></figure>
where $\rho$ is density, $u$ the velocity vector, $p$ the pressure, $\eta$ the dynamic viscosity and $F$ the volume force vector (gravity). There is no time dependency in the model as they are solved to a steady state, so the terms $\rho \frac{\partial u}{\partial t}$ and $\rho \frac{\partial p}{\partial t}$ in Eqs. (1) and (2) are neglected.

## 2.2 Magma Composition

The properties of the magma are modelled as the averaged properties of its constituents: melt, crystals and gas. For the reference model, the general composition of the melt is taken as rhyolitic, using the groundmass analysis of Montserrat dome rocks undertaken by Barclay et al. (1998). However, several melt compositions, were ultimately considered to assess the effect of melt composition on the modelled eruption dynamics (Table 2). Crystal content ($\chi_c$) and density ($\rho_c$) are fixed as we assume a constant temperature and that the conduit ascent times are orders of magnitudes faster than the time required for crystal growth by decompression, meaning only the phenocrysts present in the magma chamber are accounted for and growth of microlites and microphenocrysts is not considered. The expression for the bulk density of the magma is given by:

$$ \rho = \rho_m\chi_m(1 - \chi_g) + \rho_g\chi_g + \chi_c\rho_c(1 - \chi_g), \quad (3) $$

where $\chi_m$ is the initial fraction of melt ($1 - \chi_c$), $\rho_m$ is the melt density and $\chi_g$ is the gas volume fraction (Table 1). For the gas phase, water is assumed as the only volatile species present and the gas density ($\rho_g$) is calculated directly from the ideal gas law with the assumption that bubble growth is in equilibrium with the conduit pressure:

$$ pV = nRT, \quad (4) $$

where $V$ is the volume of gas, $R$ the ideal gas constant and $T$ the temperature. The number of moles of water, $n$, is related to density by:

$$ n = \frac{M}{m}, \quad (5) $$

where $M$ is the molar mass of water and $m$ is the mass of water present. Thus, combining Eqs. (7) and (8) and considering a unit volume we get:

$$ \rho_g = \frac{mp}{RT} \quad (6) $$

In the reference model, a single magma temperature is used with the exception of the temperature across a thermal boundary layer (TBL) defined adjacent to the conduit wall. A linear temperature drop is applied across the TBL, to simulate the cooling of the magma abutting the country rock in a well-established conduit (Collier and Neuberg 2006). The gas volume fraction ($\chi_g$) is calculated by determining how much water remains dissolved within the melt at a particular pressure using the solubility of H₂O in rhyolitic melts presented by Liu et al. (2005). At significant pressures, all the water is dissolved within the melt fraction and $\chi_g$ is initially zero. However, as the pressure decreases, water begins to exsolve out of the melt and forms bubbles. The absolute volume of exolved gas (V) can be calculated from the ideal gas law (4). This absolute volume of gas is then used to calculate the gas volume fraction of the bulk magma constituted by the gas phase.

## 2.3 Magma Viscosity—The Contribution of Crystals, Bubbles and Melt

The bulk magma viscosity ($\eta$) is determined by first calculating the viscosity of the pure melt phase ($\eta_m$). This is done using a model for the viscosity of magmatic liquids presented by Giordano et al. (2008), that predicts the viscosity of silicate melts as a function of temperature and melt composition. It is important to note that the composition used in the viscosity model is that of the pure melt phase (rhyolitic) not the overall
<table><caption>**Table 2** Compositions of melt used in the numerical simulations</caption><thead><tr><th>Composition<sup>a</sup></th><th>SiO<sub>2</sub></th><th>Al<sub>2</sub>O<sub>3</sub></th><th>TiO<sub>2</sub></th><th>FeO</th><th>MgO</th><th>MnO</th><th>CaO</th><th>Na<sub>2</sub>O</th><th>K<sub>2</sub>O</th></tr></thead><tbody><tr><td>a</td><td>71.41</td><td>13.58</td><td>0.28</td><td>2.78</td><td>1.64</td><td>0.13</td><td>4.86</td><td>3.73</td><td>1.6</td></tr><tr><td>b</td><td>76.97</td><td>11.21</td><td>0.29</td><td>1.89</td><td>0.26</td><td>0.12</td><td>1.29</td><td>4.07</td><td>2.37</td></tr><tr><td>c</td><td>77.10</td><td>9.83</td><td>0.18</td><td>1.17</td><td>0.22</td><td>0.10</td><td>1.52</td><td>4.14</td><td>1.72</td></tr><tr><td>d</td><td>78.66</td><td>11.20</td><td>0.39</td><td>1.93</td><td>0.30</td><td>0.10</td><td>1.48</td><td>3.57</td><td>2.38</td></tr></tbody></table>

<sup>a</sup>Compositions determined through, (a) rastered electron microprobe analysis of groundmass (Barclay et al. 1998); (b) Matrix glass composition (Rutherford and Devine 2003); (c) Quartz hosted melt inclusion (Devine et al. 1998); (d) Cameca SX50 microprobe analysis of interstitial glass (Burgisser et al. 2010). All melts are rhyolitic and composition (a) is used in the defined reference model

magma composition. The whole rock composition of recent Soufrière Hills Volcano magma is andesitic (Edmonds et al. 2010), but this includes the contribution of the crystals. The viscosity model is only used to calculate the actual viscosity of the liquid component (the melt), on which the crystals (the solid) have no bearing. When the effect of crystals within the melt is considered, the effective viscosity of the melt (liquid) and crystal (solid) mixture ($\eta_{mc}$) increases, and can be represented by the Einstein-Roscoe equation:

$$ \eta_{mc} = \eta_m \left( 1 - \frac{\chi_c}{\chi_c^{\max}} \right)^{-2.5}, \quad (7) $$

where $\chi_c^{\max}$ is the volume fraction of crystals at which the maximum packing is achieved and a commonly adopted value for this is 0.6 (Marsh, 1981), which is used within this study. Although this value was proposed for randomly packed spheres, and it has been shown by Marti et al. (2005) that $\chi_c^{\max}$ tends to decrease as the particle (crystal) shape becomes less isotropic. Ishibashi (2009) demonstrated that this value is a good approximation as the effect of particle shape on $\chi_c^{\max}$ is offset by effects of size heterogeneity and crystal alignment.

The presence of bubbles also affects the viscosity. If the bubbles within the magma remain un-deformed they act to increase viscosity, whilst if deformed (elongated in the direction of flow), they act to decrease viscosity (Llewellin and

Manga, 2005). Whether a bubble is in an un-deformed or deformed state is represented by the capillary number:

$$ Ca = \frac{\eta_m r E}{\Gamma} \quad (8) $$

where r is the un-deformed bubble radius, $\Gamma$, the bubble surface tension and E, a function of the strain rate within the magma flow defined below. If $Ca > 1$ then the bubbles can be considered deformed. The value of Ca will vary as a function of shear strain rate and elongation strain rate (Thomas and Neuberg 2012), meaning bubbles can be deformed within the model through either shear or extension. To account for strain acceleration or deceleration the dynamic capillary number (Cd) is required (Llewellin and Manga 2005). This compares the timescale over which the bubbles can respond to changes in their strain environment with the timescale over which the strain environment changes. If this value is large, the flow is termed unsteady and the bubbles are unable to deform independently in response to the flow. However, for the models considered here, conditions of unsteady flow are found only in a very small area near the exit of the conduit. Accounting for this within the models resulted in no noticeable change in the derived flow parameters, hence the computation of Cd is not considered.

Depending on the value of Ca, $\eta$ is calculated using the suggested 'minimum variation' of Llewellin and Manga (2005):
$$ Ca = \begin{cases} <1 & \eta = \eta_{mc}(1 - \chi_g)^{-1} \\ >1 & \eta = \eta_{mc}(1 - \chi_g)^{5/3} \end{cases} \quad (9) $$

By assuming the homogeneous nucleation of a number of bubbles in a unit volume of melt as a single event, which is determined experimentally through the initial bubble number density ($b_{ni}$) (e.g. Hurwitz and Navon 1994), the bubble radius (Lensky et al. 2002) is given by:

$$ r = \left[ \frac{S_0^3 \rho_m (C_0 - C_m)}{\rho_g} \right]^{1/3}, \quad (10) $$

where $C_0$ and $C_m$ are the initial and remaining amount of water dissolved in the melt respectively and $S_0$ is the initial size of the melt shell from which each bubble grows. $S_0$ is related to the instantaneous bubble number density ($b_n$) through the expression:

$$ S_0^3 = \frac{3}{4\pi b_n}. \quad (11) $$

$b_n$ is used rather than the initial value ($b_{ni}$) because homogeneous nucleation is assumed. Therefore, the bubble number density must remain constant with respect to the volume of the melt fraction in Eq. 3. This also accounts for bubble coalescence and $b_n$ is given by:

$$ b_n = \frac{b_{ni}}{\chi_m} [\chi_m - (1 - \chi_g)] \quad (12) $$

## 2.4 Brittle Failure of Melt

It is now well established that magma, or more specifically the melt component of a magma can fail in a brittle manner (e.g. Goto 1999). This is likely to generate low-frequency (LF) earthquakes (e.g. Neuberg et al. 2006) and effect the overall flow dynamics. In order to account for these effects, it is necessary to define conditions under which the melt may fracture. Shear failure of melt occurs when the shear stress ($\eta\dot{\epsilon}$) exceeds the

shear strength ($\tau_s$), and has been represented as a brittle failure criterion (e.g. Tuffen et al. 2003):

$$ \frac{\eta\dot{\epsilon}}{\tau_s} > 1 \quad (13) $$

where $\dot{\epsilon}$ is the shear strain rate. This criterion holds true under the assumption that during un-relaxed deformation the accumulation of shear stress in the melt obeys the Maxwell model:

$$ \sigma_s = \frac{\eta}{\mu} \frac{\partial \sigma_s}{\partial t} = \eta \dot{\epsilon} \quad (14) $$

where $\sigma_s$ is the shear stress and $\mu$ the shear modulus.

The magma composition as discussed in Sects. 2.2 and 2.3 is considered without the effects of microlite growth, the reasons for doing so are outlined in Sect. 2.2. However, it is worth noting, that if considered, the influence of microlite growth would possibly increase the bulk viscosity significantly.

## 2.5 Boundary Conditions

Flow within the system is driven by a pressure gradient defined by boundary conditions at the top and bottom of the conduit. The top boundary is set to atmospheric pressure at the altitude of the conduit exit plus any overburden load from an emplaced lava dome. The bottom boundary is set to lithostatic pressure (assuming a homogeneous country rock density of 2600 kg m⁻³) plus any imposed overpressure ($P_e$). Both the top and bottom pressure conditions are held constant throughout the model run. Initial boundary conditions along the length of the conduit are defined as no slip. When brittle failure of melt is considered within a model run, at the regions of the conduit wall where the brittle failure criterion was exceeded, the boundary conditions are changed to a tangential slip velocity ($\Delta u$) defined by:

$$ \Delta u = \frac{1}{\beta} \sigma_s, \quad (15) $$
where $\sigma_s$ is the tangential shear stress to the conduit wall and the coefficient $\beta$ is a function of the slip length ($L_s$):

$$ \beta = \frac{\eta}{L_s}, \qquad (16) $$

The model is then re-run to account for the effect of changing boundary conditions at the conduit walls. Where this results in an increase in the predicted failure depth, an iterative approach is used and the model is re-run with the new depth until the depth at which brittle failure of the melt stabilises. For the purposes of this study, the failure depth is considered to have converged if the depth increase between iterative runs is less than 10% of the previous observed increase.

The size and shape of the conduit is also an important factor in influencing the accent rate, so for the purpose of assessing the potential magnitude of this influence, several possible conduit shapes were modelled. While the types of boundary conditions discussed above do not change, the relative locations of the boundaries do (Fig. 2). Case (a) is the simplest geometry change and represents just a change of the conduit radius (r). Case (b) represents a narrowing of the conduit. Case (c) represents a widening of the conduit. The extent to which the geometry of the conduit is changed within the models is discussed further in the next Section.

# 3 Critical Conduit Processes

## 3.1 Using Magma Ascent Rates to Assess Model Sensitivity

The ascent rate is a key parameter in understanding volcanic hazard because it has been directly linked to eruptive behaviour (e.g. Gonnerman and Manga 2007). By gaining a better understanding of which model parameters have the greatest effect on ascent rates, we can achieve an insight into which are the most important parameters controlling explosivity, and the likely severity of the volcanic hazard. For the purpose of comparing the various models we use two velocities, defined as $\bar{V}$ and $V_{2500}$, where $\bar{V}$ is the average accent velocity taken along a vertical profile through the centre of the conduit, and $V_{2500}$ is the average accent velocity taken along a horizontal profile at a depth of 2500 m within the conduit.

The ascent rate has also been linked to monitoring data such as seismicity (e.g. Thomas and Neuberg 2012) or deformation (e.g. Zobin et al. 2011), therefore, it is possible to link the changes in model parameters to recorded monitoring data. In addition, there are physically observed variations in ascent rate estimated from a variety of methods, ranging from studying mineral reaction rims around phenocrysts within erupted magma

<figure><img src="image_59.png" /><figcaption>**Fig. 2** Schematic diagram of the alternative conduit geometries modelled, showing a a constant conduit radius; b a narrowing conduit and c a widening conduit</figcaption></figure>
(Rutherford and Devine 2003) to interpreting lava dome morphology (Sparks et al. 2000). This places constraints on the magnitude of changes to the modelled ascent rate engendered by altering the model input parameters we can consider realistic.

Matching the absolute values of physically observed and calculated ascent rates is currently beyond the scope of the model, however we can use the magnitude of the observed variations to provide upper and lower bounds to the extent to which the model input parameters are varied. Any changes that produce increases in ascent rate greater than two orders of magnitude over the reference model are not considered realistic in this work. This may seem at first an arbitrary discrimination, but there is a good reason that the observed or calculated ascent rates presented in the literature (e.g. Rutherford and Devine 2003; Castro and Gardner 2008) are “slow” ($<5 \times 10^{-2} \text{ m s}^{-1}$). Faster ascent rates, while likely to exist in nature, would almost certainly result in substantial fragmentation of the magma, making it very difficult to observe or calculate the actual magma ascent rate below the initial point of fragmentation. Fragmentation dynamics are not considered within the current model, hence no valid inferences or conclusions can be gained from studying the model runs that exhibit extremely fast ascent rates.

## 3.2 The Critical Model Parameters

Figure 3 summarises the sensitivity of ascent rate to the different model parameters presented in Table 1. The single parameter (within the modelled ranges) which has the strongest effect on the ascent velocities is the initial dissolved water content of the magma. This parameter affected both $\bar{V}$ and $V_{2500}$ to a large degree. In contrast there are several model parameters which have little effect on the modelled ascent velocities. These include the thermal boundary layer thickness and the temperature drop across it, as well as the bubble number density and bubble surface

tension. Modifying the parameters involved in the brittle failure of the melt (magma shear strength and slip length) has a negligible effect on ascent rates and these results have not been plotted on Fig. 3. However, the contribution of the brittle failure of the melt to observed geo-physical signals is considered very important, and will be discussed in Sect. 3.3.

It is unsurprising that the group of model parameters that appear to have the greatest effect on the magma ascent velocity, as seen in Fig. 3b (water content, temperature, crystal content, and chemical composition), also have the greatest effect on the magma viscosity (Sect. 2.3). Ultimately, modelling the ascent of magma is a fluid flow problem, and the properties that have the biggest effect on the fluid (magma) properties will have the biggest effect on the overall dynamics of the system. All other parameters have a much smaller direct effect on the fluid properties, and although they may be important to specific small scale magmatic processes when considered in isolation, with respect to the overall magma ascent they appear insignificant. For example, altering the properties of the bubbles within the magma, $b_{ni}$ and $\Gamma$, the effect is to change the shape and the number of bubbles. Previous work has heavily focused on this area (e.g. Llewellin and Magna 2005) but the effect on the overall flow modelled here is minimal. The indication from this is that it is the total volatile content (water in this case) which is available that is more important to governing the overall flow dynamics, rather than how exactly it is stored in the magma. This particular observation is a key point as new estimates from Cassidy et al. (2015) suggest that basaltic South Soufrière Hills magmas (and by extension, possibly other basaltic arc magmas) have the potential to be extremely volatile-rich, containing up to >6 wt% H$_2$O prior to eruption. Firstly, this validates the use of high initial water contents used in the range of parameters modelled, and secondly, given the range of ascent velocities generated within the models as a result of just changing the dissolved water content (the dark blue bars in
<figure><img src="image_60.png" /><figcaption>**Fig. 3 a** The ascent rate, $\bar{V}$ for each value of the parameters altered. The black bar represents the reference model with parameters as listed in Table 1. **b** The same data as (a) plotted relative to the reference model, which is represented by the horizontal black line. Parameter</figcaption></figure>

changes that caused increases in the ascent rates plot upward from the horizontal line, while parameter changes that caused decreases in ascent rates plot downwards. **c** $V_{2500}$ for each value of the parameters altered, plotted relative to the reference model

Fig. 3), demonstrates that it is vitally important to obtain an accurate understanding of the magma components at the volcano of interest, rather than assuming “typical” values representative of a broad compositional category.

Outside of water content, temperature, crystal content, and chemical composition, the two

parameters modelled which have the largest influence on the modelled ascent rate are the chamber overpressure and the conduit geometry. These are particularly important points when considering volcanoes entering periods of unrest following long periods of quiescence. It is problematic to achieve an accurate understanding
of the magma components highlighted above at all volcanoes under these circumstances due to a likely lack of monitoring (a problem highlighted in Parts 1 and 2 of this book). Unless there has been long-term measurement of deformation occurring at the volcano now exhibiting signs unrest, it will be extremely difficult to estimate any likely overpressures in the chamber, and attempting to define the conduct geometry of a system that has not yet erupted would be almost impossible. It is therefore paramount that as much information as possible of all potentially active volcanic systems is routinely gathered before signs of unrest are detected.

## 3.3 Matching Observations—Explosivity and Seismicity

Although, as previously mentioned, matching the absolute values of physically observed and calculated ascent rates is currently beyond the scope of the model, key to giving the models real significance is determining whether the changes to important parameters highlighted in Sect. 3.1 can be theoretically linked to physical observations at real volcanic systems. Figure 4 shows values of $V_{2500}$ for all of the modelled parameters in a manner similar to that presented in Fig. 3c, but in this case, the data are plotted relative to a baseline accent value of $0.02 \text{ ms}^{-1}$. This baseline

value was chosen because it has been highlighted by Rutherford and Devine (2003) as an ascent rate which may indicate a transition between effusive and explosive behaviour. This value has been obtained from quantifying the breakdown of hornblende in ascending magma, and while this technique is not an accurate barometer for defining an exact ascent velocity required for explosive eruptions, the rates calculated for non-explosive eruptive activity at Soufrière Hills volcano between the period of November 1995–September 2002 were below this value. Figure 4 shows that several model runs produced ascent rates of $<0.02 \text{ ms}^{-1}$ (by altering the melt composition or $\chi_c$) and several other runs produced ascent rates very close to this value (by altering $W_{\epsilon c}$, T, and conduit geometry), indicating that by altering just single parameters within the system this theoretical threshold of ascent rate can be crossed.

Conduit flow is treated as a closed system, so no outgassing is considered in the model, as a result the ascent velocities are overestimated (Thomas and Neuberg 2014). It is therefore predictable that if this process was included, far more of the model runs would result in ascent velocities that straddle the baseline in Fig. 4. This suggests that the ability for the ascent rate within the conduit to fluctuate, either side of values that have been linked to explosive eruptions in response to small changes in the system

<figure><img src="image_61.png" /><figcaption>**Fig. 4** The ascent rate, $V_{2500}$ for each value of the parameters altered, plotted relative to an ascent velocity of $0.02 \text{ ms}^{-1}$</figcaption></figure>
parameters is genuine. The requirement to accu- rately understand and model outgassing pro- cesses (see Sect. 4) is therefore an important capacity that is currently lacking.

One major discrepancy between physical observations and the model is the model results suggested that the brittle failure of the melt (Sect. 2.4) and the related LF seismicity would occur as a shallow process (in agreement with the work of Holland et al. (2011)). The physical observations place the location of this type of seismicity at Soufrière Hills at depths of ~1500 m below the conduit exit (Neuberg et al. 2006). This is because under normal conditions the shear stress required to break the melt (Eq. 14) can only be reached where the melt is extremely viscous, which occurs near the surface. In order to reach the higher shear stresses required to break the melt at greater depths, where the viscosity is lower, the shear strain rate (\dot{\epsilon}) needs to increase. Since \dot{\epsilon} is equal to the lateral velocity gradient within a cylindrical conduit or dyke:

$$ \dot{\epsilon} = \frac{dv}{dx} \qquad (17) $$

the simplest way to increase \dot{\epsilon} is to increase the velocity of the magma flowing within the conduit, or reduce the area through which it flows,

which since mass must be conserved also has the effect of increasing the flow velocity.

To resolve this discrepancy between model and observations we introduce a constriction within the conduit as a plausible explanation for brittle failure at greater depths. We test its effect within the reference model by including a bottleneck region at a depth of 1500 m, reducing the conduit diameter from 15 to 10 m. This bottleneck is 100 m in length, which equates to only 1/50 of the total conduit length. Figure 5 shows ascent velocity and shear strain rate profiles from the bottleneck region compared to values from the same location of the conduit in the unmodified reference model. By altering this relatively small region of the conduit, the shear strain rate increases by a factor of four. Crucially, with the exception of small changes in the magma rheology caused by the induced pressure gradients within the bottleneck, the magma viscosity has not been altered. Due to the increased value of shear strain rate the brittle failure ratio (13) will increase by the same factor. By introducing such asperities into the conduit and increasing the strain rate it is possible to drive the brittle fracture of the melt to deeper levels in the conduit that match the location of recorded LF seismicity at Soufrière Hills volcano. This further

<figure><img src="image_62.png" /><figcaption>**Fig. 5** a Plot of the shear strain rate within a simple bottleneck of 100 m length at intervals of 0.005. The values of shear strain rate are seen to be increases over the entire length of the bottleneck. b Cross conduit profiles</figcaption></figure>

<figure><img src="image_63.png" /><figcaption>taken at the same depth for velocity and shear strain rate for the case of the unmodified reference model (solid line) and a conduit containing a bottleneck (dashed line)</figcaption></figure>
emphasises the importance of understanding the possible conduit geometry previously highlighted in Sect. 3.2. This is just one possible solution of generating seismicity deeper in the conduit. Any process that acts to increase the sear strain rate, such as strain localisation between crystals in a magma with a high crystal fraction (as discussed in Chap. 10 of this book) would be a possible explanation.

## 4 Pathways for Outgassing

The role of volatiles has been identified as the primary controlling parameter governing ascent rate within the conduit, but what was not considered was that the gas, once exsolved from the melt has the ability move independently to the melt. This is a commonly assumed process in basaltic systems, but not in silicic magmas of a relatively high viscosity as it is assumed that bubbles of gas cannot rise with significant speed within the magma. However fractures generated by brittle failure of the melt (Sect. 3.3) may provide ideal outgassing pathways for exsolved gas.

The behaviour of this exsolved gas is separated from the problem of magma ascent, and considered independent of any other parameters, through additional numerical modelling in Comsol Multiphysics (Collinson and Neuberg, 2012). In order to consider all possible outgassing pathways, including vertically through the conduit, or laterally through the walls, we model the gas response to brittle failure using a simplistic “block-style” model, with a central conduit and adjacent wall-rocks. Brittle failure is explicitly modelled as an increase in the permeability within narrow regions either side of the conduit. The problem is simplified to considering permeable flow through a static media. Consequently, the equations for Darcy’s law (18), and the continuity equation (19) are amalgamated to derive a partial differential equation (20), which is solved for pressure (P):

$$ \mathbf{u} = -\frac{k}{\eta_g}(\nabla P + \rho g \nabla z) \quad (18) $$

$$ \frac{\partial}{\partial t}(\rho\varepsilon) + \nabla \cdot (\rho\mathbf{u}) = 0 \quad (19) $$

$$ \frac{\partial}{\partial t}(\rho\varepsilon) + \nabla \cdot \rho \left[ -\frac{k}{\eta_g}(\nabla P + \rho g \nabla z) \right] = 0 \quad (20) $$

The gas velocity (u) is then determined by using the pressure gradient within Darcy’s law (18).

In contrast, with the models for ascent rate, time dependency is included in this model in order to understand the changes in the system through time, in response to a permeability (k) increase in response to brittle failure at the conduit-wall margin. Due to its abundance in volcanic systems, water is the only volatile considered here. The gas density (ρ) is calculated using the mean molar mass within the ideal gas law (6), and the gas viscosity (μ) assumed constant at $1.5 \times 10^{-5}$ Pa (Collinson and Neuberg 2012).

Bulk permeabilities are set such that the conduit has a higher permeability ($10^{-10}$ m²) than the wallrocks ($10^{-12}$ m²), and initially, the conduit-wall margin is “sealed” with a very low permeability of $10^{-16}$ m². The fracturing is initiated at 1500 m, in accordance with measurements by Neuberg et al. (2006), and propagate vertically towards the surface, as an increase in permeability from $10^{-16}$ to $10^{-6}$ m².

Figure 6a shows the initial system, before fracturing, where the gas loss is predominantly vertical, through the conduit. In Fig. 6b, the fracture zone has propagated vertically up to 700 m depth. Consequently, the pressure has increased through the conduit and corresponding wall margins where the fractures have developed. This is due to the regions of increased permeability being confined by areas of lower permeability above. Thereby, providing a suitable environment for gas storage, which due to the increased pressurisation, may force exsolved volatiles back into solution within the melt. This change in pressurisation has resulted in a
<figure><img src="image_64.png" /><figcaption>**Fig. 6** **a** The initial system in equilibrium, before fracturing has commenced. **b** The systems of fractures has propagated upwards to a depth of 700 m. **c** The fracture zones have reached the surface. **d** The system has</figcaption></figure>

resumed equilibrium. Velocity range of arrows (ms⁻¹): **a** max: 0.30, min: 1.3 × 10⁻³; **b** max: 0.38, min: 8.8 × 10⁻⁵; **c** max: 0.33, min: 6.8 × 10⁻⁴; **d** max: 0.30, min: 9.5 × 10⁻⁴

corresponding change to the gas velocity pattern, which the lower conduit having a very low gas velocity due to the low pressure gradient. However, at shallow levels, the higher pressure gradient forces a much higher gas velocity towards the surface. On reaching the surface (Fig. 6c), the high permeability fractures result in a decrease in pressure throughout the conduit and a rapid expulsion of the stored gas. Equilibrium conditions resume approximately 4 h after the fracture zones reached the surface (Fig. 6d). This model shows a propagating fracture zone is an effective mechanism for degassing the conduit and wall margins. A key observation is the possibility for this mechanism to produce periods of cyclic activity which are observed at many silicic volcanoes (e.g. Holland et al. (2011)), which can be directly related to observed degassing patterns or through controlling the ascent rate and the

associated geophysical signals (e.g. LF seismicity) within the conduit through moderating the amount of gas stored in the system.

# 5 Summary and Implications

In the introduction to this chapter, we asked the question of whether it is possible that small changes in the composition of properties of the magma could cause changes significant enough to be recorded in monitoring data or even visual observations. It is clear that these small changes, particularly when considering changes in the water content or conduit diameter, can have large effects of the ascent velocity of the magma. These effects are large enough that conceivable they could be simply observed as an increased exclusion rate at the surface. At this point however it
may be too late. If increased extrusions rates are being observed at the surface a critical threshold of ascent rate may have already been surpassed.

More useful would be to observe these potential changes through monitoring before the magma was physically observed to be extruding faster at the surface. Throughout this chapter we have discussed the importance of shear stress, seismicity and outgassing, and it is to monitoring data relating to these processes that we would look for an indication of potential changes.

Increases in ascent velocity may lead to a change in eruption style, but they also potentially cause changes in monitoring data. As discussed in the chapter, any increase in ascent velocity of the magma will cause an increase in the shear stress experience by the magma, and will likely lead to an increase in rate of seismicity, or alter its location. It is now generally accepted that shear stress within the conduit also causes a significant deformation signal (e.g. Neuberg et al. 2018). Any changes in the number or location of low frequency earthquakes, or changes in the near-field deformation around the conduit could therefore be inferred as changes in the magma properties. In addition, the single most important parameter identified within this study was changes in the water content. The fractures generated by the seismicity, as a result of the increases in ascent velocity, have also been shown to be important outgassing pathways. An increase in water content would lead to an increase in volatiles and these would be more easily outgassed along the created fractures. Conceptually, a clear signal of an increase in ascent velocity (and a potential early warning of a change in eruption style) caused by increase in water content of the magma, could be an increase and deepening of low frequency earthquakes, accompanied by an increase in the deformation signal followed by significant outgassing event.

In this chapter we look at a single volcanic setting, but by changing the model parameter these effects could be assessed at any volcano, and regardless of the volcano studied the relative importance of the parameters considered should remain unaltered.

**Acknowledgements** The authors are grateful for inspiring discussions with Geoff Kilgour (GNS New Zealand) during J. Neuberg's study leave at GNS Research Centre in Wairakei, and contributions to research ideas form former Ph.D. students within the Volcanic Studies Group at the University of Leeds.

## References

Barclay J, Rutherford M, Carroll M, Murphy M, Devine J (1998) Experimental phase equilibria constraints on pre-eruptive storage conditions of the Soufrière Hills magma. Geophys Res Lett 25:3437–3440

Burgisser A, Poussineau S, Arbaret L, Druitt TH, Gischetti T, Bourdier J-L (2010) Pre-explosive conduit conditions of the 1997 Vulcanian explosions at Soufrière Hills Volcano, Montserrat: I. Pressure and vesicularity distributions. J Volcanol Geoth Res 194:27–41

Castro JM, Gardner JE (2008) Did magma ascent rate control the explosive-effusive transition at the Inyo volcanic chain, California? Geology 36:279–282

Collier L, Neuberg J (2006) Incorporating seismic observations into 2D conduit flow modelling. J Volcanol Geoth Res 152:331–346

Collinson ASD, Neuberg JW (2012) Gas storage, transport and pressure changes in an evolving permeable volcanic edifice. J Volcanol Geoth Res 243–244:1–13

Cluzel N, Laporte D, Provost A, Kannewischer I (2008) Kinetics of heterogeneous bubble nucleation in rhyolitic melts: implications for the number density of bubbles in volcanic conduits and for pumice textures. Contrib Miner Petrol 156:745–763

Devine JD, Murphy MD, Rutherford MJ, Barclay J, Sparks RSJ, Carroll MR, Young SR, Gardner JE (1998) Petrologic evidence for pre-eruptive pressure-temperature conditions, and recent reheating, of andesitic magma erupting at the Soufrière Hills Volcano, Montserrat, W.I. Geophys Res Lett 25:3669–3672

Devine JD, Rutherford MJ, Norton GE, Young SR (2003) Magma storage region processes inferred from geochemistry of Fe–Ti oxides in andesitic magma, Soufrière Hills Volcano, Montserrat, WI. J Petrol 44:1375–1400

Edmonds M, Aiuppa A, Humphreys M, Moretti R, Giudice G, Martin RS, Herd RA, Christopher T (2010) Excess volatiles supplied by mingling of mafic magma at an andesite arc volcano. Geochem Geophys Geosyst 11:Q04005

Giordano D, Russell JK, Dingwell DB (2008) Viscosity of magmatic liquids: a model. Earth Planet Sci Lett 271:123–134

Gonnermann HM, Manga M (2007) The fluid mechanics inside a volcano. Annu Rev Fluid Mech 39:321–356
Gotto A (1999) A new model for volcanic earthquake at Unzen Volcano: melt rupture model. Geophys Res Lett 26:2541–2544

Holland ASP, Watson IM, Phillips JC, Caricchi L, Dalton MP (2011) Degassing processes during lava dome growth: Insights from Santiaguito lava dome, Guatemala. J Volcanol Geoth Res 202:153–166

Hurwitz S, Navon O (1994) Bubble nucleation in rhyolitic melts: experiments at high pressure, temperature, and water content. Earth Planet Sci Lett 122:267–280

Ishibashi H (2009) Non-Newtonian behavior of plagioclase-bearing basaltic magma: subliquidus, viscosity measurement of the 1707 basalt of Fuji volcano, Japan. J Volcanol Geoth Res 181:78–88

Lensky NG, Lyakhovsky V, Navon O (2002) Expansion dynamics of volatile-supersaturated liquids and bulk viscosity of bubbly magmas. J Fluid Mech 460:39–56

Liu Y, Zhang Y, Behrens H (2005) Solubility of H₂O in rhyolitic melts at low pressures and a new empirical model for mixed H₂O-CO₂ solubility in rhyolitic melts. J Volcanol Geoth Res 143:19–235

Llewellin E, Manga M (2005) Bubble suspension rheology and implications for conduit flow. J Volcanol Geoth Res 143:205–217

Lyakhovsky V, Hurwitz S, Navon O (1996) Bubble growth in rhyolitic melts: experimental and numerical investigation. Bull Volc 58:19–32

Marsh B (1981) On the crystallinity, probability of occurrence, and rheology of lava and magma: contributions to mineralogy and petrology. Bull Volcanol 78:85–98

Marti I, Höfler O, Fischer P, Windhab EJ (2005) Rheology of concentrated suspensions containing mixtures of spheres and fibers. Rheol Acta 44:502–512

Neuberg JW, Tuffen H, Collier L, Green D, Powell T, Dingwell D (2006) The trigger mechanism of

low-frequency earthquakes on Montserrat. J Volcanol Geoth Res 153:37–50

Neuberg JW, Collinson ASD, Mothes PA, Ruiz CM, Aguaiza S (2018) Understanding cyclic seismicity and ground deformation patterns at volcanoes: Intriguing lessons from Tungurahua volcano. Ecuador, Earth Planet Sci Lett 482:193–200

Rutherford MJ, Devine JD (2003) Magmatic conditions and magma ascent as indicated by hornblende phase equilibria and reactions in the 1995–2002 Soufrière hills magma. J Petrol 44:1433–1454

Sparks RSJ, Murphy MD, Lejeune AM, Watts RB, Barclay J, Young SR (2000) Control on the emplacement of the andesite lava dome of the Soufrière Hills volcano, Montserrat by degassing-induced crystallization. Terra Nova 12:14–20

Thomas ME, Neuberg JW (2012) What makes a volcano tick—A first explanation of deep multiple seismic sources in ascending magma. Geology 40:351–354

Thomas ME, Neuberg JW (2014) Understanding which parameters control shallow ascent of silicic effusive magma. Geochem Geophys Geosyst 15:4481–4506

Tuffen H, Dingwell DB, Pinkerton H (2003) Repeated fracture and healing of silicic magma generate flow banding and earthquakes? Geology 31:1089–1092

Wadge G, Robertson REA, Voight B (eds) (2014) The eruption of Soufrière Hills Volcano, Montserrat from 2000 to 2010. Geol Soc, London, Memoirs, p 39

Zobin VM, Ramirez JJ, Santiago H, Alatorre E, Navarro C (2011) Relationship between tilt changes and effusive-explosive episodes at an andesitic volcano: the 2004–2005 eruption at Volcan de Colima, Mexico. Bull Volc 73:91–99

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

<figure><img src="image_65.png" /></figure>
<figure><img src="image_66.png" /></figure>

# When Does Magma Break?

Fabian B. Wadsworth, Taylor Witcher, Jérémie Vasseur,
Donald B. Dingwell and Bettina Scheu

## Abstract

Geophysical signals arriving at the Earth's surface originate from a source mechanism at depth but are not necessarily directly observable. Therefore, well-posed experiments can provide insights into source mechanics and, importantly, the parameters required to model aspects of the sources of unrest signals. In this Chapter we detail one such example of how experimental laboratory work has improved our understanding of unrest signals. We focus on the failure of single- and multi-phase magmas, demonstrating that the liquid viscosity, and therefore the temperature and volatile content of a magma of a given composition, is the limiting parameter in determining whether a magma will ascend viscously or whether it can fracture during ascent. This critical threshold is characterized by a Deborah number, the ratio of the timescale of relaxation to the timescale of local flow. We show that for single-phase magmatic liquids and for vigorously vesiculating magmas, a local Deborah number of $10^{-2}$ is the limit above which mixed viscoelastic behaviour including fracture propagation can be expected, and a Deborah number of 1 is the limit above which magma is dominantly elastic and responds in a brittle manner to applied stresses. These thresholds can be understood in terms of the onset and peak of the Debye relaxation process for viscoelastic liquids. The apparent validity of a Maxwell model permits us to predict the maximum stress that can be supported by a volcanic liquid deforming in the high Deborah number range. We use these constraints to provide a map of timescales on which we contour dominant system responses from viscous to purely brittle; valid for all magmatic liquids. Finally, we explore the scaling necessary to extend these conceptual insights to crystal- and bubble-bearing magmas valid under specific conditions.

<footnote>F.B. Wadsworth (✉) · T. Witcher · J. Vasseur
D.B. Dingwell · B. Scheu
Ludwig-Maximilians-Universität, Theresienstr.41,
80333 Munich, Germany
e-mail: fabian.wadsworth@min.uni-muenchen.de</footnote>
The competing timescales of deformation and relaxation in magma are relevant to unrest source mechanisms that originate from magma deformation, such as long-period seismic signals that are used to predict eruption timing.

**Keywords**

Rheology · Strain rate · Experimental volcanology · Glass transition · Failure forecasting · Low frequency earthquakes · Viscous dissipation

**Glossary**

<table><tr><td>Rheology</td><td>The study of the response of a material to an applied stress or deformation. In the volcano-sciences this typically refers to magma-rheology which is dominated by the rheology of the liquid phase (a viscous or viscoelastic melt) and the additional effect of suspended phases (crystals and gases).</td></tr><tr><td>Viscoelasticity</td><td>A material response to an applied stress in which both elastic and viscous components of the deformation are observed. In magma, as the temperature decreases or the local strain rate increases, the elastic component becomes more dominant. When the viscous component is negligible, purely elastic behaviour and material-rupture can readily occur. See Rheology.</td></tr></table>

# 1 Introduction

There are a wide variety of observable unrest signals at active volcanoes. The key unrest signals are those that are diagnostic of a new regime of behaviour or of how the system may evolve in the future. One such family of events are the low-frequency earthquakes at volcanoes, which are thought to directly represent magma movement at moderate to shallow depths (Chouet et al. 1994; Chouet 1996; Neuberg et al. 2006; Salvage and Neuberg 2016). Such events are powerful tools particularly because they are closely associated with surface activity (Miller et al. 1998). Accelerating event rates of low frequency earthquakes can be used to retrospectively

predict an eruption or dome collapse event with reasonable accuracy (Salvage and Neuberg 2016). The utility of the low-frequency event type has been solidified since physical hypotheses for the source mechanism have been put forward (Neuberg et al. 2006) and tested in scaled laboratory experiments (Benson et al. 2008; Tuffen et al. 2008). The phenomenological observation of the use of certain signal types over others is useful and has led to accurate retrospective eruption forecasts (Voight 1988; Chouet 1996). However, knowledge of the physics of the source leads to a more diverse range of tools being deployed to understand the evolution of volcanic systems into the future. In the context of low-frequency earthquakes at volcanoes, once a
source mechanism is identified (Tuffen et al. 2003; Neuberg et al. 2006), numerical models of magma ascent in the conduit beneath volcanoes can be used to reproduce the depths of the source from first principles (Thomas and Neuberg 2012) and could eventually make forecasts of other observables at the surface that would be consistent with impending eruption.

In this Chapter we explore the source mechanism of low-frequency earthquakes at volcanoes from a physical perspective using a compilation of data from scaled laboratory experiments. We use these datasets to demonstrate that the critical threshold for fracture propagation in a viscoelastic fluid such as magma is universal. We propose that this potentially simplifies the inputs required for effective modelling of source mechanics of low-frequency events. More than this, we hope that this provides a good example of how laboratory work can provide valuable insight into the physical feasibility of models proposed to explain unrest signals at volcanoes.

# 2 Scaling the Viscous-to-Brittle Transition in Magmas

## 2.1 Single-Phase Magmatic Liquids

Newtonian viscous liquids deform and flow under applied stresses such that the rate of shear strain $\dot{\gamma}$ resulting from an applied shear stress $\tau$ is proportional to the viscosity $\mu$ via $\tau = \mu\dot{\gamma}$. However, as large shear stresses are applied to viscous liquids, they can exhibit viscoelastic behaviour, which manifests as an apparent non-Newtonian relationship between $\tau$ and $\dot{\gamma}$ and can result in fracture propagation locally in the liquid. Of interest here is the transition from simple, apparently Newtonian viscous behaviour to complex viscoelastic behaviour.

Maxwell described the simplest viscoelastic model in which a liquid has a characteristic time required to relax an applied shear stress, termed the Maxwell relaxation timescale $\lambda_r$, which is given by $\lambda_r = \mu/G_\infty$ where $G_\infty$ is the shear modulus in the purely elastic regime. The Maxwell model of viscoelasticity permits us to

predict that if a shear strain is imposed on a viscoelastic liquid, the resultant shear stress will rise rapidly to a peak value $\tau_i$ and will relax over time $t$ according to $\tau(t) = \tau_i \exp(-t/\lambda_r)$. This simple concept has proved powerful in describing the transition between Newtonian and viscoelastic behaviour in volcanic liquids (Dingwell 1995, 1996). Indeed, this model is invoked to explain the transition between volcanic liquids and volcanic glasses on cooling (e.g. Stevenson et al. 1995), the fragmentation of fluid and bubbly magma undergoing decompression (Alidibirov and Dingwell 1996; Kameda and Kuribara 2008) and the conditions under which extensive shear fracture networks can form at or near volcanic conduit margins (Gonnermann and Manga 2003; Tuffen et al. 2003; Kendrick et al. 2014; Hornby et al. 2015). In all cases, it is useful to define a Deborah number De, which is the dimensionless ratio between the Maxwell relaxation time and the timescale of deformation $\lambda$. For simple (viscometric) shearing flow, the latter timescale is $1/\dot{\gamma}$ and thus

$$De = \frac{\lambda_r}{\lambda} = \frac{\mu\dot{\gamma}}{G_\infty} \quad (1)$$

Our hypothesis is that a Deborah number of unity separates the regime between viscous (De $\ll$ 1) and elastic (De $\gg$ 1) responses to shear stresses, or equivalently, between a system that can relax shear stresses efficiently from one in which shear stresses can accumulate elastically. In two different experimental types, both Webb and Dingwell (1990a) and Cordonnier et al. (2012a, b) showed that in fact the first evidence of behaviour that is not purely viscous occurs at De $\ge$ $10^{-2}$ and not exactly at unity. This observation will be discussed later.

The power of Eq. 1 is that simple parameters that can be measured in the laboratory can be used to predict where this transitional point of De = 1 will be. This critical Deborah number can be termed De'. In what follows, we give examples of how $\mu$ and $G_\infty$ have been found for a range of volcanic liquid compositions at magmatic temperatures and volatile contents before exploring the validity and applicability of the
scaling in Eq. 1. The dimensionless nature of Eq. 1 means that it can be assessed for any system so long as $\lambda$ can be defined.

At the core of the experimental toolkit is the determination of the fundamental quantities that are required to model volcanic processes involved in unrest. The most variable physical quantity in volcanic systems is the viscosity of the liquid phase. Using a range of techniques, this viscosity can be determined with prodigious accuracy and significant effort has been expended in mapping the full range of composition, temperature (see Giordano et al. 2008) and volatile content (Hess and Dingwell 1996) relevant to shallow magmatic settings. Multi-component models have been proposed such that the viscosity of any composition of silicate volcanic liquid on Earth can now be predicted as a function of temperature including the range of conditions between shallow magma storage and the surface (Hess and Dingwell 1996; Giordano et al. 2008). This provides parameterization of $\mu(T)$ for use in Eq. 1. The most commonly used form for $\mu(T)$ is the non-Arrhenian Vogel-Fulcher-Tammann equation of the form

$$ \mu = A \exp\left(\frac{B}{T-C}\right), \quad (2) $$

where A, B, and C are coefficients that are experimentally determined and then parameterized as a function of families of oxides in the liquid structure (Giordano et al. 2008) or as a function of the dissolved water content (for calc-alkaline rhyolites: Hess and Dingwell 1996), to give two examples.

To give examples of the range of viscosities of interest in volcanic scenarios, we give end members in Fig. 1 for a calc-alkaline rhyolite using the model from Hess and Dingwell (1996), and for the basaltic liquid composition provided in Zhang et al. (1991) using the model from Giordano et al. (2008), both contoured for a range of water contents.

Unlike viscosity, the shear modulus in the elastic regime does not vary significantly in silicate liquids and is not strongly dependent on composition or temperature. Indeed, a short survey of the values for $G_{\infty}$ in silicate glasses and liquids at a range of temperature and a huge range of composition, provided by Dingwell and Webb (1989) shows that $G_{\infty} = 10^{10 \pm 0.5}$ Pa. Here we use this range in order to fully parameterize De (by predicting $\lambda_r$) as a function of temperature (via Eq. 2) for any silicate liquid in the shallow crust.

## 2.2 Extensions to Multiphase Magmas

Except in rare circumstances such as obsidian-forming eruptions, volcanic liquids are not often erupted without some proportion of suspended pore space (either as isolated bubbles or connected networks) and rigid crystals. In either case, the utility of Eq. 1 requires additional attention. We posit that for the liquid phase between the pores or the crystals, De given by Eq. 1 holds. However, we acknowledge that estimation or measurement of the rate of shear strain locally between pores or crystals would be difficult (Deubelbeiss et al. 2011). Therefore, a more robust criterion for the viscous-to-brittle transition in crystal- or bubble-bearing magma would require that we scale the bulk strain rate $\dot{\gamma}_b$ on the system for the effect of the suspended load.

For crystals, we show some first-order scaling attempts for their effect on the critical threshold for the onset of brittle behaviour in multiphase magma. The simplest view of the local flow of liquid in crystal-bearing magma under constant shear stress is that the rate of shear strain between the crystals should scale approximately with $\dot{\gamma} = \dot{\gamma}_b(1 - \phi_x/\phi_m)^{-1}$ where $\phi_x$ is the suspended crystal volume fraction and $\phi_m$ is a jamming fraction above which no more crystals can be
<figure><img src="image_67.png" /><figcaption>**Fig. 1** The viscosity of end-member magmatic liquids with variable dissolved water concentrations. Plotted are the results for calc-alkaline rhyolite from the general viscosity model for hydrous silicic liquids (Hess and Dingwell 1996) for 0.01–1 wt% water, and a typical basaltic composition (composition from Zhang et al.</figcaption></figure>

1991) calculated using the multicomponent viscosity model (Giordano et al. 2008) for 0–1 wt% water. The relaxation timescale is calculated assuming a composition independent value for $G_\infty$ of $10^{10}$ Pa (Dingwell and Webb 1989)

added to the flowing system. This scaling would imply that the Deborah number for a crystal-bearing magma $De_x$ would be

$$De_x = \frac{\mu}{G_\infty} \dot{\gamma}_b \left(1 - \frac{\phi_x}{\phi_m}\right)^{-1}, \quad (3)$$

where $\phi_m$ is a function of crystal shape and roughness (Mueller et al. 2010; Mader et al. 2013). A critical value of $De_x$, termed $De'_x$, is $10^{-2}$, consistent with the limiting $De'$ at $\phi = 0$. A new definition of the bulk failure criterion $De'$ can now be made, which decreases as $\phi \to \phi_m$ so that $De' = De'_x(1 - \phi_x/\phi_m)$. In Fig. 2 we demonstrate how this concept predicts a linear relationship between $De_x$ and $\phi_x/\phi_m$, which in turn shows that lower bulk strain rates are required to induce brittle behaviour in crystal-bearing magmas.

This hypothesis was found to hold for simple two-phase systems by Cordonnier et al. (2012a), but where those authors required that $\mu$ be replaced by the suspension viscosity of the whole

system, which is at odds with the scaling of strain rate for local liquid effects only. Nonetheless, this approach described their data satisfactorily and the conceptual insight that crystals locally increase the liquid rate of shear strain relative to the bulk value is robust, with the implication is that the whole suspension will begin fracturing at lower bulk rates of strain relative to a single-phase liquid of the same composition as crystals are added. It may be that additional second-order effects are important at high $\phi_x$, which are not accounted for in this simple analysis.

For the case of porous magmas, there are two considerations: (1) The growth of bubbles exerts a rate of shear strain in the liquid concentrated at the bubble walls, which is broadly independent of any bulk shearing deformation, and (2) like for the crystal case, the presence of bubbles changes the partitioning of the bulk rate of shear strain in the liquid between the bubbles. In the case of scenario (1), the Rayleigh-Plesset equation can be used to relate the bubble growth rate to the gas
<figure><img src="image_68.png" /><figcaption>**Fig. 2** The critical Deborah number required for fracturing $De'$ is reduced when crystals are suspended in magma. To a first-order, the $De'$ value is reduced proportional to $De' = De'_x(1 - \phi_x/\phi_m)$</figcaption></figure>

pressure in the bubble, relative to the hydrostatic pressure (Sparks 1978). This can be augmented for the diffusion-controlled gradient of viscosity in the immediate liquid shell around a growing bubble (Prousevitch et al. 1993; Lensky et al. 2001). The component of the rate of shear strain tangential to the bubble wall $\dot{\gamma}_{\theta}$ can then be computed throughout bubble growth (Ichihara et al. 2002). If we were then to input $\dot{\gamma}_{\theta}$ into our computation of $De$, we would predict under which conditions the system would meet the criterion of $De' = 10^{-2}$ locally at the bubble rim. These conditions would be best cast in terms of the critical rate of pressure or temperature change, or bulk initial volatile content that would allow the bubble to grow sufficiently fast that fractures could be propagated at the bubble wall.

In case (2), the effect of bubbles on the bulk Deborah number is less clear. We speculate that at $De \gg 1$, the stress required to fracture a bubbly system will be analogous to that required to fracture a vesicular glass. In this case, models for the effect of spherical cavities on the stress required for fracturing are valid and they predict that the stress is reduced significantly as the bulk gas volume fraction increases (Sammis and Ashby 1986). This has been confirmed in the high Deborah number regime for porous liquids analogous to volcanic systems (Vasseur et al.

2013) but remains untested in the low Deborah number regime where bubble deformation may be important.

## 2.3 Apparent Non-newtonian Effects

In single-phase liquids at Deborah numbers below the critical threshold at which fracturing is observed, there is evidence that during steady state shearing flow, there is a non-Newtonian relationship between applied shear stress and resultant shear rate of strain (Simmons et al. 1982; Dingwell and Webb 1989; Webb and Dingwell 1990a; Cordonnier et al. 2012b). In the onset of this non-Newtonian, the onset of this non-Newtonian behaviour appears to be well described by the value above which viscous dissipation of heat is active on the system length scale of interest (Costa and Macedonio 2005; Hess et al. 2008). This can be scaled by the Brinkman number $Br$

$$Br = \frac{\Phi_g}{\Phi_l} = \frac{\mu\dot{\gamma}^2}{kq}, \quad (4)$$

where $\Phi_g$ and $\Phi_l$ are the gain and loss power densities, respectively, $q$ is the areal heat flux out of the system and $k$ is the thermal conductivity of
the material. $\Phi_g = \mu\dot{y}^2$ represents the amount of energy produced in a system of a given volume by viscous dissipation of heat, while $\Phi_l = kq$ represents the energy lost due to diffusive thermal equilibration. When $Br \gg 1$, heat is efficiently produced and inefficiently lost from the system, resulting in a bulk temperature increase in the liquid. This would manifest itself as an apparent shear thinning rheology if the temperature increase where not locally accounted for, and would be most likely to be operative at high viscosities and high shear strain rates. We note that unlike the Deborah number, the Brinkman number is scale dependent (as $q$ depends on the area available for heat transfer out of the system) and so should be assessed for each system scale separately.

## 3 The Universal Breaking Timescales of Volcanic Liquids

Laboratory data related to the viscous-to-brittle transition in magmas has been collected in a variety of geometries. Using single-phase liquids, there are two dominant geometries: (1) thin fibers of silicate liquid of basaltic, andesitic, phonolitic and rhyolitic composition were stretched under constant load in tension until the fibers snapped in a singular fracture event (Webb and Dingwell 1990b), and (2) cylinders of synthetic borosilicate liquid were compressed under constant load and the bulk temporal evolution of the rate of axial shortening was determined as viscous (relaxed) if it continuously increased to a steady or near-steady value, or brittle (unrelaxed) if the axial shortening rate jumped due to fracturing events (Cordonnier et al. 2012b). In another type of experiment, analogue vesicular liquids were decompressed at different rates from pressure (Kameda and Kuribara 2008; Kameda et al. 2013). In this type of decompression experiment, if the dominant response was viscous (relaxed), then the sample was seen to grow due to decompression-driven bubble growth, and if the dominant response was brittle (unrelaxed), then

smooth sample inflation was punctuated by visible fractures opening. In these decompression experiments there additionally was a violent rupture mode in which the fracturing was pervasive and shattered the sample in a vigorous fragmentation event. These data are selected here because the liquids used have a known relaxation time under the conditions used in the experiments and care was taken by the authors who originated the work to investigate the shear rates of strain local to the bubbles (discussed above).

In Fig. 3 we map these experimental results as a ratio of the deformation timescale to the relaxation timescale, which permits us to contour the plot for critical Deborah numbers. We colour-code the data according to the bulk mode of response of the sample to the deformation using green to represent purely viscous relaxed behaviour, orange to represent brittle unrelaxed behaviour and red to represent the complete violent rupture of the sample. These data suggest that the viscous relaxed behaviour transitions to unrelaxed brittle behaviour at a Deborah number of $10^{-2}$, and are consistent across a huge range of experimental conditions and across both the single-phase compression and vesicular decompression experiment types. Furthermore, we see that a Deborah number of unity consistently separates the experiments for which the bulk response was unrelaxed and brittle from those for which the response was violent rupture, fragmentation or complete failure. It appears that these two thresholds, $De = 10^{-2}$ and $De = 1$ are universal even when comparing analogue room temperature liquids with high temperature silicate liquids deformed in a variety of ways. This lends power to the scaling provided by the dimensionless Deborah number and implies that we need only to define the liquid viscosity and a characteristic rate of deformation in order to predict whether a system will flow viscously or rupture violently. For example, the working viscosity of the analogue fluid used in Kameda and Kuribara (2008) is $10^0 \le \mu \le 10^{10}$ Pa.s, which extends to much lower used in Cordonnier et al. (2012b), and yet the scaling with the Deborah number
<figure><img src="image_69.png" /><figcaption>**Fig. 3** A map of behaviour arising from deformation experiments on natural (Webb and Dingwell 1990b) and synthetic (Cordonnier et al. 2012b) silicate liquids and on analogues of magmatic liquids undergoing vesiculation (Kameda and Kuribara 2008; Kameda et al. 2013). The experimental results were obtained by applying a characteristic timescale of deformation $\lambda$ of a liquid with a characteristic timescale of stress relaxation $\lambda_r$. Marked are ratios $\lambda_r/\lambda$, equivalent to Deborah numbers $De$, of $10^{-2}$</figcaption></figure>

and 1, which separate experiments with a purely viscous response from those with an unrelaxed brittle response or a violent rupture response, respectively. Additionally marked is the Brinkman number $Br$ of 1, which separates the boundary between isothermal viscous behaviour (above the line) and viscous behaviour in which the material heats up due to viscous dissipation as heat (below the line).

holds across a broad range simply because the deformation timescale was also smaller in the former example.

In the viscous field, Cordonnier et al. (2012b) additionally recorded whether samples hosted a measurable temperature increase due to viscous dissipation of heat during the experiment, or not. In Fig. 3, we plot those with a measureable temperature increase as filled symbols and those without this feature as unfilled symbols. To explain this, we plot the threshold dimensionless Brinkman number of unity $Br = 1$, using an estimation of the loss power density for the furnace and sample size used $\Phi_l = 10^{4.5} \text{W.m}^{-2}$ (Cordonnier et al. 2012b). We note that the Brinkman number curve consistently divides the regimes of purely isothermal experiments and those with measurable heat gain. On this map, the position of $Br = 1$ is non-unique and depends on sample size, such that on the scale of a volcanic conduit, for example, $Br = 1$ would occur

at much higher deformation timescales for the same relaxation timescale (the curve would shift up in Fig. 3). This implies that in the natural case, viscous dissipation of heat may be far more important than shown here for the sample length scale (Costa and Macedonio 2005; Mastin 2005; Costa et al. 2007). Nevertheless, the Deborah number limits discussed appear to be universal and, importantly, are scale-independent.

To gain helpful physical insight into why $10^{-2} \le De \le 10^0$ is the transitional window between a purely viscous and a purely elastic response of a liquid to a deformation, we provide low-strain, high frequency oscillatory rheological measurement data for similar liquids (Fig. 4). Here, rods of single-phase liquid are subjected to a low-amplitude oscillatory strain with a forcing frequency $\omega$ at a range of temperatures similar to those used in experiments presented in Fig. 2. Here, $\omega$ is normalized with $\lambda_r$, yielding a dimensionless frequency or, equivalently, a