# Storage and Eruption of Silicic Magma across the Transition from Dominantly Effusive to Caldera-forming States at an Arc Volcano (Santorini, Greece)

Gareth Nicholas Fabbro, Timothy H. Druitt, Fidel Costa

► To cite this version:

Gareth Nicholas Fabbro, Timothy H. Druitt, Fidel Costa. Storage and Eruption of Silicic Magma across the Transition from Dominantly Effusive to Caldera-forming States at an Arc Volcano (Santorini, Greece). Journal of Petrology, 2017, 58 (12), pp.2429 - 2464. 10.1093/petrology/egy013. hal-01781381

HAL Id: hal-01781381

https://uca.hal.science/hal-01781381v1

Submitted on 25 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Storage and Eruption of Silicic Magma across the Transition from Dominantly Effusive to Caldera-forming States at an Arc Volcano (Santorini, Greece)

Gareth N. Fabbro<sup>1,2*</sup>, Timothy H. Druitt<sup>1</sup> and Fidel Costa<sup>2</sup>

<sup>1</sup>Laboratoire Magmas et Volcans, Université Clermont Auvergne, CNRS, IRD, OPGC, F-63000 Clermont-Ferrand, France; <sup>2</sup>Earth Observatory of Singapore, Nanyang Technological University, N2-01B-30, 50 Nanyang Drive, Singapore 639798, Singapore

*Corresponding author. Present address: Earth Observatory of Singapore, Nanyang Technological University, N2-01B-30, 50 Nanyang Drive, Singapore 639798, Singapore. Telephone: +65 6592 1803. E-mail: gfabbro@ntu.edu.sg

Received February 4, 2016; Accepted February 12, 2018

## ABSTRACT

Caldera-forming eruptions at Santorini discharge large volumes of silicic magma from upper crustal reservoirs. Sequences of smaller eruptions preceding the main explosive eruptions can provide insight into the conditions of the plumbing system that lead to caldera-forming events, which is important for interpreting monitoring data. We analysed textures, zoning patterns, trace element compositions, and crystal residence timescales calculated at pre-eruptive conditions (near-liquidus residence timescales; NLRT) of plagioclase and orthopyroxene phenocrysts from two eruptive units on Santorini: an ~2 km³ sequence of dacitic lavas erupted between 39 and 25 ka (the Therasia dome complex), and a caldera-forming dacitic eruption that occurred no more than 2800 ± 1400 years after the last Therasia lava (the 21.8 ka, >10 km³ Cape Riva eruption). The study builds on our previous work, which showed that the Therasia and Cape Riva dacites, although similar in most major elements, differ in some trace element contents and were derived from different source magmas and crystal mushes. Contents of K and La in plagioclase phenocrysts mirror those in the respective host magmas (higher in Therasia, lower in Cape Riva), showing that plagioclase provenance can be determined using these elements. Very few plagioclase cores from the Cape Riva dacite were recycled from Therasia magmas; the majority were derived from lower-K, lower-La magmas and mushes related to the Cape Riva eruption itself. Despite the very different magma volumes and eruptive fluxes, plagioclase and orthopyroxene crystals from the two dacite series have remarkably similar textures and major element compositions. Furthermore, Mg diffusion profiles in plagioclase and Mg-Fe diffusion profiles in orthopyroxene yield similar ranges of NLRT, most ranging from years to centuries. Some orthopyroxene crystals exhibit Al sector zoning indicative of rapid growth. Processes driving crystallization appear to have been similar in the two systems, despite the differences of scale. Based on previously published phase diagrams and melt inclusion volatile barometry for the Cape Riva dacite, we infer that in each case crystallization of plagioclase rims and orthopyroxene took place centuries to years prior to eruption owing to volatile-saturated decompression (± cooling) as the dacitic melts (plus entrained plagioclase antecrystic cores) ascended from the middle crust (10–16 km) into the upper crust (4–6 km), where they resided until eruption a few years to decades later. Between 39 and 25 ka, multiple small volumes of Therasia-type dacitic magma were emplaced in the upper crust, where they either froze or were subsequently erupted. From about 25 ka onwards, large volumes of Cape Riva-type dacitic magma, sourced from a different mid-crustal reservoir, began to ascend into the upper crust. Runaway
drainage of this magma source, peaking during the decades to years prior to the Cape Riva erup-
tion, led to establishment of a well-mixed magma chamber in the upper crust that was discharged
during the caldera-forming Cape Riva event.

**Key words:** calderas; crystal zoning; interplinian–Plinian transition; magma reservoirs; Santorini

INTRODUCTION

Caldera-forming eruptions of silicic magma are a major hazard to society and the environment (e.g. Self & Blake, 2008). They are commonly preceded by long periods of relatively low eruptive flux involving effusive activity, weak explosive activity and dormancy (e.g. Metz & Mahood, 1985; Stix et al., 1988). It is during such periods that the conditions are established in the crust for the next caldera-forming eruption. Better understanding of those conditions, and of the processes involved, may help us anticipate the outcome of future unrest at caldera volcanoes. In this paper we use information recorded in crystals from a prehistoric caldera-forming eruption at Santorini, and from a sequence of lavas and tuffs that preceded it, to reconstruct the assembly, then explosive discharge, of a caldera-forming silicic magma reservoir.

The magma reservoirs that feed caldera-forming eruptions are believed to be located in the upper crust (e.g. Scaillet & Evans, 1999; Liu et al., 2005; Hildreth & Wilson, 2007; Gertisser et al., 2009; Arce et al., 2013; Cadoux et al., 2014). Silicic melts at Santorini are probably generated in deep crustal hot zones by fractional crystallization and partial melting, then transferred to the shallow crustal storage reservoir (Annen et al., 2006; Solano et al., 2012; Allan et al., 2013, 2017; Cadoux et al., 2014; Andújar et al., 2015, 2016; Druitt et al., 2016). At many volcanoes such as Santorini, caldera-forming eruptions involve magmas that have few crystals (typically <20 vol. %). In these cases, the melt must have been separated from its source crystal mush and segregated prior to ascent and eruption (Hildreth, 1981; Hildreth & Wilson, 2007; Bachmann & Bergantz, 2008).

One model for the construction of large upper crustal reservoirs is that they are assembled incrementally over long periods. Comparing radiometric ages for systems with high closure temperatures (such as U-Pb in zircon or Rb-Sr in feldspar) with estimates of eruption age has provided evidence for long crystal residence times, suggesting that some melts can be stored for up to 10⁴–10⁵ years prior to their eruption (e.g. Christensen & DePaolo, 1993; Davies et al., 1994; Brown & Fletcher, 1999; Vazquez & Reid, 2004; Folkes et al., 2011; Wotzlaw et al., 2013). Effusive eruptions of compositionally similar magma prior to caldera-forming eruptions have been interpreted as recording reservoir growth over periods of 10⁴–10⁵ years (Bacon, 1985). Proposed examples of ‘precursory leaks’ from such reservoirs include the Glass Mountain rhyolites that were erupted

prior to the 760 ka Bishop Tuff eruption (Hildreth, 1979; Metz & Mahood, 1985, 1991; Halliday et al., 1989; Hildreth & Wilson, 2007), the pre-climactic lavas that were erupted prior to the 6850 years BP climactic eruption of Mount Mazama, Oregon (Bacon & Druitt, 1988; Druitt & Bacon, 1989; Bacon & Lowenstern, 2005), and the Tihoi and Okaia eruptions prior to the 26.5 ka, 530 km³ Oruanui eruption at Taupo, New Zealand (Sutton et al., 1995).

Recent work suggests, however, that reservoir assembly may occur much more rapidly. Chemical gradients in phenocrysts found in the eruptive products of some large silicic eruptions have been interpreted as documenting late-stage recharge events, and diffusion models indicate that these gradients cannot have existed for more than 10²⁻¹⁰³ years at near-liquidus temperatures (Wark et al., 2007; Boyce & Hervig, 2008; Saunders et al., 2010; Druitt et al., 2012; Gualda et al., 2012; Wilcock et al., 2012; Allan et al., 2013). Other work has cast doubt on whether ‘precursory leaks’ tap the same reservoir as the subsequent caldera-forming eruptions. Zircon populations in the Glass Mountain rhyolites are not the same as those in the Bishop Tuff (Simon & Reid, 2005; Chamberlain et al., 2014). Similarly, different zircon populations found in the Okaia, Tihoi and Oruanui eruptions of Taupo show that the three melt bodies were assembled separately in the shallow crust shortly before their respective eruptions, although they may have shared a deeper reservoir over a longer period of time (Wilson & Charlier, 2009). The Oruanui magma can only have been emplaced in the upper crust after the Okaia eruption, which gives a residence time in the shallow reservoir of less than 3000 years, as also supported by the study of Mg-Fe diffusion gradients preserved in orthopyroxene crystals in the Oruanui pumice (Allan et al., 2013, 2017).

The 21.8 ± 0.4 ka, >10 km³ Cape Riva caldera-forming silicic eruption of Santorini, Greece was preceded, from 39 to 25 ka, by a period of relatively low eruptive flux (<0.1 km³ ka⁻¹) during which an ~2 km³ sequence of silicic lavas and minor tuffs up to 200 m thick called the Therasia dome complex was emplaced (Fabbro et al., 2013). Despite similarities in major element bulk-rock chemistry, the Therasia dacites differ in their trace element signatures from the Cape Riva dacite, and were probably derived from a different parent (Fabbro et al., 2013). High-resolution ⁴⁰Ar/³⁹Ar dating shows that the last Therasia lava was extruded no earlier than 24.6 ± 1.3 ka, so the sequence offers a record of the evolution of the crustal magmatic system over most of the 18 kyr preceding the Cape Riva eruption. In this
paper, we use major and trace element zoning in plagio-
clase and orthopyroxene phenocrysts from the
Therasia and Cape Riva products to investigate the
processes that took place during the 18 kyr build-up to
the Cape Riva eruption. We use the zoning patterns of
slow-diffusing elements to calculate melt compositions
in equilibrium with the crystals, and the gradients of
fast-diffusing elements to constrain crystal residence
timescales. The new data are combined with the previ-
ous constraints of Fabbro et al. (2013) to reconstruct the
evolution of the Santorini magmatic system across the
transition from dominantly effusive to caldera-forming,
explosive states.

VOLCANOLOGICAL BACKGROUND

Santorini lies on the Hellenic volcanic arc, formed by
the subduction of the African plate beneath the Aegean
Sea (Le Pichon & Angelier, 1979; Jolivet et al., 2013), on
continental crust about 25 km thick (Tirel et al., 2004;
Karagianni et al., 2005; Endrun et al., 2008). Volcanic ac-
tivity started on Santorini at ~650 ka, and since 360 ka
there have been at least 12 Plinian eruptions leading to
at least four caldera collapses (Druitt et al., 1999).
Interplinian periods were characterized by effusive and
minor explosive activity, and multiple lava edifices have
been constructed at different times inside, or around
the edge of, the caldera: the Peristeria cone (530-
430 ka), the Skaros shield (67–54 ka), the Therasia dome
complex (39–25 ka), and the present-day Kameni edifice
(<3.6 ka; Druitt et al., 1999).

Santorini magmas range from basalt to rhyodacite,
are transitional tholeiitic to calc-alkaline in character,
and since 530 ka almost all have had anhydrous pheno-
cryst assemblages (Huijsmans et al., 1988; Huijsmans &
Barton, 1989; Druitt et al., 1999). Most Santorini melts
are saturated with an H₂O-rich fluid at upper crustal
pressures (Cadoux et al., 2014; Druitt et al., 2016).
Intermediate to silicic melts form principally by frac-
tional crystallization of basaltic parents (50 wt % SiO₂,
7 wt % MgO) containing 1–4 wt % H₂O in the melt
(Andújar et al., 2016; Druitt et al., 2016). Phase equilibria
studies of a basalt (Andújar et al., 2015), an andesite
(Andújar et al., 2016), a dacite, and three rhyodacites
(Cadoux et al., 2014), along with H₂O + CO₂ melt inclu-
sion barometry (Druitt et al., 2016), serve to constrain
the architecture of the plumbing system beneath the
caldera and its evolution over time. Basalt differentiates
to andesite at about 16 km depth (400 MPa), close to the
limit between the upper and lower crust (Sachpazi et al.,
1997; Karagianni et al., 2005; Endrun et al., 2008;
Konstantinou, 2010), and andesite fractionates to silicic
compositions at about 10–16 km depth (250–400 MPa).
Andesitic and silicic melts ascend into the upper crust,
where they are stored at depths of a few kilometres
prior to final eruption (Cadoux et al., 2014; Druitt et al.,
2016).

The focus of this study is the 21.8 ± 0.4 ka Cape Riva
eruption, and the effusive eruptions of dacitic lava that

preceded it (Druitt, 1985; Fabbro et al., 2013; Figs 1 and
2). Between 39.4 ± 2.2 ka and 24.6 ± 1.3 ka, extrusion of
at least 12 dacitic lava domes and flows (hereafter sim-
ply referred to as flows), interspersed with minor explo-
sive eruptions, formed the ~2 km³ Therasia dome
complex (Fabbro et al., 2013). The vents of the dome
complex were located across the western flank and
summit of a basaltic-to-andesitic lava shield, called the
Skaros shield, which had grown between ~70 and 54 ka
(Huijsmans, 1985; Huijsmans & Barton, 1989; Druitt
et al., 1999). The summit of the Skaros-Therasia edifice
subsequently collapsed during the Cape Riva eruption
(Druitt, 1985; Druitt et al., 1999; Fig. 1).

The dacites of the Therasia lavas and Cape Riva
pumice both have 64–68 wt % SiO₂ (Fig. 2), and contain
up to 20 vol. % crystals of plagioclase, orthopyroxene,
clinopyroxene, magnetite, ilmenite and apatite, and

<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Morphological evolution of Santorini between 54 and 21 ka, after Druitt et al. (1999). The dashed line is the present-day outline of the islands. Contours are at 100 m intervals. (a) The maximum extent of the Skaros shield, constructed between ~70 and 54 ka. (b) The maximum extent of Therasia dome complex at ~25 ka. (c) The island shortly after the ~22 ka Cape Riva eruption.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Summary of magma compositions and ages between ~70 and 20 ka, after Fabbro *et al.* (2013). Major explosive eruptions are represented by stars, with periods of edifice construction coloured in yellow. Periods of little or no preserved eruptive activity are left blank. The SiO₂ contents of the eruptive products are shown on the right. Black symbols are whole-rock compositions, and white symbols are groundmass compositions. For the Cape Riva eruption, the large stars represent the dacite that makes up >95% of the deposit, and the small star is a hybrid andesite that is found mingled with the dacite in the Plinian unit (Cape Riva A) of that eruption. Brown bars represent paleosols. The Therasia and Cape Riva chemical data are taken from Fabbro *et al.* (2013), and the Skaros chemical data from Huijsmans (1985). ¹Dates from Fabbro *et al.* (2013); ²dates from Schwarz (2000); ³dates from Druitt *et al.* (1999).</figcaption></figure>

trace amounts of xenocrystic olivine (Druitt, 1985; Fabbro *et al.*, 2013). Fe–Ti oxide thermometry gives 880 ± 25°C for both dacites (Fabbro *et al.*, 2013; Cadoux *et al.*, 2014). The Therasia dacites contain quenched enclaves of basaltic composition, as well as dispersed fragments of sub-solidus olivine-bearing gabbro. The last erupted lava prior to the Cape Riva eruption was an olivine-bearing hybrid andesite ('Upper Therasia andesite'; Fig. 2) formed by the mixing of approximately equal proportions of Therasia dacite and basalt. Dacite makes up >99% of the magmatic components of the Cape Riva eruption, but small blebs of quenched olivine-bearing basalt occur in pumices from the first eruption phase, as well as a hybrid andesitic magma formed by mixing of the dacite and basalt that is found as individual clasts and in banded pumices mingled with the dacite.

Despite similarities in most major element contents, the Cape Riva dacite (hereafter called CR dacite) is poorer in incompatible elements [K, La, Rb, Zr; also Nb, Ta, Th, Hf and light rare earth elements (LREE)] than the Therasia dacites (hereafter called TH dacites) at a given SiO₂ content (Fabbro *et al.*, 2013). The same difference is observed between the respective hybrid andesitic

components, showing that the two dacites also had different basaltic parents (referred to as TH and CR basalts). The decrease in incompatible elements from the Therasia series to Cape Riva series is one step in a long-term temporal trend towards lower incompatible elements at Santorini over the last 530 kyr. Although the decrease is most evident in silicic magmas, it is also evident in the respective intermediate and basaltic members of each series (Huijsmans *et al.*, 1988; Bailey *et al.*, 2009; Fabbro *et al.*, 2013). The chemical differences between different basaltic parents at Santorini is attributed to different degrees of partial melting and of metasomatism of the mantle source by subducted sediment (Huijsmans *et al.*, 1988; Bailey *et al.*, 2009; Vaggelli *et al.*, 2009). Magma evolution within each series can be modelled by fractional crystallization with minor crustal assimilation (Bailey *et al.*, 2009).

The model proposed by Fabbro *et al.* (2013) involves the following stages: (1) the construction of the high Skaros lava shield from ~70 to 54 ka, which culminated in the eruption of the large, explosive, andesitic Upper Scoria 2 eruption; (2) a pause between 54 and 40 ka, during which only a few andesitic lavas were erupted and most basalt ascending from the mantle was trapped at depth; (3) the slow leakage of ~2 km³ of incompatible-enriched TH dacite (with sparse basaltic enclaves and hybrid andesites) from 39 to 25 ka; (4) the first possible minor appearance of incompatible-poorer CR dacite in the upper crust at 26 ka; (5) the accelerated arrival of CR dacite (with basalt) into the shallow crust between 25 and 22 ka, followed by its caldera-forming eruption at 21.8 ± 0.4 ka. Assembly of the upper crustal Cape Riva magma reservoir took place within the final 2800 ± 1400 years prior to its eruption.

## ANALYTICAL TECHNIQUES

Plagioclase and orthopyroxene crystals were studied in three Therasia dacitic lavas [flows 5, 8 and 9 of Fabbro *et al.* (2013)] and in six Cape Riva pumice lumps [two from phase A, one from phase B, two from phase C and one from phase D of Druitt (1985)]. The three Therasia lavas were selected from the section with the best stratigraphic control, to span as much of the period of effusive, dacitic activity as possible. The Cape Riva samples were chosen to cover all four phases of the eruption. Crystals in the lavas were studied in polished thin sections. Those in the pumices were extracted by crushing and set in polished grain mounts; crystals with adhering melt were favoured, as this shows that the surfaces are original crystal faces.

Mineral major element compositions for >180 crystals were analysed using the Cameca SX 100 electron microprobe at the Laboratoire Magmas et Volcans, Université Clermont-Auvergne, Clermont-Ferrand, using a beam current of 15 nA, a voltage of 15 kV and a focused beam (nominally 1 µm in diameter).

Backscattered electron (BSE) images of >300 plagioclase crystals were produced on the scanning electron
microscope (SEM), and a subset of these were cali-
brated for anorthite content [An = 100 × molar Ca/
(Na + Ca + K)] using the electron microprobe. Five
points on each crystal were analysed for An content
using the electron microprobe, and the greyscale values
of these points were measured on the SEM images
using the image analysis software package ImageJ
(Rasband, 1997). Correlation between greyscale value
and An content had coefficients *r*² > 0.97 for most of the
images (and *r*² > 0.90 for all). This corresponds to a stand-
ard error in the calculated An content of ± 1–2 mol %. The
images were smoothed using a 2 pixel radius median fil-
ter. The advantage of a median filter over a moving aver-
age is that it removes random noise while leaving any
genuine sudden variations in composition unchanged.
This means that sharp boundaries between crystal rims
and cores were not softened by the smoothing process.
We also smoothed the An profiles during the numerical
modelling of diffusion to ensure stability of the solution
(Costa *et al.*, 2003) but this has no significant impact on
our results or interpretation.

The zoning patterns of >200 orthopyroxene crystals
were examined using BSE images, and the greyscale
values of a subset of these were calibrated individually
for Mg number [Mg# = 100 × molar Mg/(Mg + Fe)].
Correlations between greyscale values and Mg# of the
orthopyroxene were generally high (r² > 0.95 for most,
and r² > 0.80 for all). This corresponds to a standard
error in Mg# of ±0.2–1 mol %. No filter was applied to
the orthopyroxene BSE images. Element maps were
made of 13 orthopyroxene crystals for Mg, Fe, Ca and
Al using a beam current of 100 nA and a dwell time of
100 ms per pixel (1–3 µm).

Trace elements in plagioclase (Li, Mg, K, Ti, Fe, Sr,
Ba, La, Ce, and Pb) were measured at the Laboratoire
Magmas et Volcans, Université Clermont-Auvergne,
Clermont-Ferrand, France using a Resonetics M-50-E
193 nm laser ablation (LA) system connected to an
Agilent 7500cs inductively coupled plasma mass spec-
trometry (ICP-MS) system. The beam diameter was 10–
15 µm and the laser pulse rate was 2 Hz. Detection limits
ranged between 0.1 and 1 ppm for most trace elements,
but were higher for K (10–40 ppm), Ti (5–20 ppm) and Fe
(30–100 ppm). Calcium was used as the internal stand-
ard, NIST610 was used as the calibration standard, and
NIST612 and BCR2G were analysed during each set of
measurements to check the calibration. For most of the
analyses, Ca was measured by electron microprobe be-
fore LA-ICP-MS analysis. Photomicrographs of the ana-
lysed crystals, with the laser pits visible, were
superimposed onto the calibrated BSE images. Relative
errors based on counting statistics for the trace element
concentrations were between 3 and 10%.

PHENOCRYST COMPOSITION AND ZONING

We use the term phenocryst for a macrocryst of which a
significant part, including the rim, formed by crystalliza-
tion from melt compositionally similar to that in which

the crystal was erupted. The term antecryst is used in
the sense of Hildreth & Wilson (2007) to mean a crystal
inherited from a genetically related magmatic precursor
(e.g. the mush from which the host melt was extracted).
The term xenocryst is used to denote a foreign crystal
or crystal component. Phenocrysts may have antecrys-
tic or xenocrystic cores.

Plagioclase

Plagioclase constitutes 75–85 vol. % of the phenocryst
content of the Therasia lavas (total crystal content 1–
17 wt %) and ~70 vol. % of that of the Cape Riva pumice
(total 15–20 wt %) (Druitt et al., 1999; Fabbro et al., 2013;
Cadoux et al., 2014). Plagioclase phenocrysts in the
Therasia and Cape Riva dacites have broadly similar
textures (Fig. 3) and compositions (Table 1). Most are
unzoned or weakly oscillatory zoned, with variations of
<5 mol % An (Fig. 3a). Two zones (core and rim) are
common to all zoned phenocrysts: (1) a compositionally
and texturally variable core, consisting either of
multiple normally zoned layers separated by resorption
surfaces ('sawtooth zoning'), or of partially resorbed,
high-An plagioclase with complex oscillatory or patchy
zoning (Fig. 3); (2) a well-defined rim zone, normally
zoned to a euhedral margin, that has the same compos-
itional range as unzoned phenocrysts in the same rock.

The core is commonly separated from the rim by a
resorption surface. In some cases, the entire zoned crys-
tal consists of a series of concentrically stacked saw-
tooth layers, the outermost of which we refer to as the
rim. In this sense, the rim records the progressive con-
vergence of the pressure, temperature and other envir-
onmental conditions towards the pre-eruptive
conditions represented by the euhedral margin of the
rim.

Plagioclase cores in the Therasia dacites range from An₉₁ to An₃₉; rims are normally zoned from An₅₉–₇₂ towards euhedral margins that vary between An₃₉ and An₅₁, depending on the lava flow. Only xenocrysts with adhering, microlite-rich groundmass derived from the disintegration of mafic enclaves have rims of An>60. There is little compositional variation of the phenocryst margins within single lava flows (<4 mol % An; Fig. 4), with the exception of occasional calcic xenocrysts, suggesting that each erupted dacite was homogeneous. However, the distribution of margin An contents varies from flow to flow; for example, rim edges in lava flow 5 are on average less calcic than those in flows 3, 8 and 9 (Fig. 4).

Zoned plagioclase grains in the Cape Riva pumice
(from all four eruptive phases) have cores ranging from
An₆₈ to An₃₆. Rims of zoned crystals (and the compos-
itions of unzoned crystals) lie within a very narrow
range of An contents (An₃₆–₄₀), with no obvious trend
with stratigraphic height (Fig. 4). Zoning profiles of
three trace elements (Sr, Mg, and La) in two representa-
tive zoned plagioclase crystals (from the Therasia and
Cape Riva dacites) are shown in Fig. 5.
<figure><img src="image_6.png" /><figcaption>Representative plagioclase phenocrysts from the Cape Riva dacite</figcaption></figure>

<figure><img src="image_7.png" /></figure>

**Fig. 3.** Calibrated backscatter electron images of plagioclase crystals from the Therasia and Cape Riva dacites, along with corresponding sketches. Anorthite spot values from electron microprobe analyses are marked on the sketches in dark blue. From the Therasia dacites: (a) a weakly zoned phenocryst (GS10-16 XL105); (b) a phenocryst without a calcic core, but with an internal resorption surface and a normally zoned rim (GS10-27a XL49); (c) a phenocryst with a glass inclusion-rich core, overgrown by a normally zoned rim (GS10-17 XL57); (d) a phenocryst with a simply zoned calcic core (GS10-14 XL66); (e) a phenocryst with a partially resorbed, complexly zoned core, overgrown by an oscillatory zoned rim (GS10-17 XL60). From the Cape Riva dacite: (f) a phenocryst with a simply zoned calcic core and an unzoned, more sodic rim (GS11-34a XL09); (g) a phenocryst with a partially resorbed calcic core, surrounded by a more sodic rim with weak oscillatory zonation (GS11-39b XL29); (h) a phenocryst without a calcic core, but with two internal resorption surfaces, each overgrown by a normally zoned sawtooth zone (GS11-34a XL14); (i) a phenocryst with complexly zoned, partially resorbed core and a normally zoned rim (GS11-30b XL08).
<table><caption>**Table 1:** Typical plagioclase compositions</caption><thead><tr><th></th><th colspan="2">Therasia (GS10-17 XL60)</th><th colspan="2">Cape Riva (GS11-34a XL14)</th></tr><tr><th></th><th>Rim</th><th>Core</th><th>Rim</th><th>Sawtooth peak</th></tr></thead><tbody><tr><td colspan="5"><strong>Electron microprobe (wt %)</strong></td></tr><tr><td>SiO<sub>2</sub></td><td>58.77</td><td>47.83</td><td>58.69</td><td>53.39</td></tr><tr><td>Al<sub>2</sub>O<sub>3</sub></td><td>25.11</td><td>31.53</td><td>25.99</td><td>29.29</td></tr><tr><td>FeO</td><td>0.45</td><td>0.54</td><td>0.36</td><td>0.47</td></tr><tr><td>MgO</td><td>0.05</td><td>0.02</td><td>0.03</td><td>0.06</td></tr><tr><td>CaO</td><td>8.16</td><td>16.05</td><td>8.16</td><td>12.17</td></tr><tr><td>Na<sub>2</sub>O</td><td>6.55</td><td>2.21</td><td>6.54</td><td>4.43</td></tr><tr><td>K<sub>2</sub>O</td><td>0.47</td><td>0.07</td><td>0.29</td><td>0.13</td></tr><tr><td>Total</td><td>99.57</td><td>98.25</td><td>100.06</td><td>99.93</td></tr><tr><td>An</td><td>79.7</td><td>39.7</td><td>40.1</td><td>59.9</td></tr><tr><td colspan="5"><strong>Cations based on 8 oxygens</strong></td></tr><tr><td>Si</td><td>2.645</td><td>2.234</td><td>2.642</td><td>2.421</td></tr><tr><td>Al</td><td>1.332</td><td>1.736</td><td>1.369</td><td>1.565</td></tr><tr><td>Fe</td><td>0.017</td><td>0.021</td><td>0.013</td><td>0.018</td></tr><tr><td>Mg</td><td>0.003</td><td>0.001</td><td>0.002</td><td>0.004</td></tr><tr><td>Ca</td><td>0.393</td><td>0.803</td><td>0.391</td><td>0.591</td></tr><tr><td>Na</td><td>0.286</td><td>0.200</td><td>0.567</td><td>0.389</td></tr><tr><td>K</td><td>0.027</td><td>0.004</td><td>0.017</td><td>0.008</td></tr><tr><td colspan="5"><strong>LA ICP-MS (ppm)</strong></td></tr><tr><td>Li</td><td>11.5</td><td>14.8</td><td>10.7</td><td>5.60</td></tr><tr><td>Mg</td><td>214</td><td>221</td><td>151</td><td>302</td></tr><tr><td>K</td><td>1830</td><td>3130</td><td>2220</td><td>904</td></tr><tr><td>Ti</td><td>212</td><td>160</td><td>—</td><td>—</td></tr><tr><td>Fe</td><td>3430</td><td>2530</td><td>2050</td><td>2600</td></tr><tr><td>Rb</td><td>0.78</td><td>1.02</td><td>—</td><td>—</td></tr><tr><td>Sr</td><td>370</td><td>383</td><td>337</td><td>352</td></tr><tr><td>Zr</td><td>2.51</td><td>0.97</td><td>—</td><td>—</td></tr><tr><td>Ba</td><td>83.5</td><td>226</td><td>168</td><td>48.3</td></tr><tr><td>La</td><td>3.29</td><td>4.95</td><td>3.80</td><td>1.44</td></tr><tr><td>Ce</td><td>5.78</td><td>6.96</td><td>5.14</td><td>2.61</td></tr><tr><td>Nd</td><td>2.52</td><td>1.82</td><td>1.47</td><td>0.58</td></tr><tr><td>Pb</td><td>24.6</td><td>4.33</td><td>2.53</td><td>0.94</td></tr></tbody></table>

Sawtooth peak refers to the most calcic part of the sawtooth zones.

## Orthopyroxene

Orthopyroxene (opx) accounts for up to ~10 vol. % of the phenocryst assemblages. Those in the Therasia dacites span the range $Wo_{3-4}En_{53-68}Fs_{29-43}$ with Mg# 58–71 [end-members calculated according to Deer et al. (1992)]. Those from the Cape Riva dacite are similar, but extend to more Fe-rich compositions ($Wo_{3-4}En_{52-65}Fs_{32-44}$; Mg# 53–68). $Al_2O_3$ and CaO contents of all opx are 0.4–1.7 and 1.4–1.9 wt %, respectively (Table 2).

Three types of opx are distinguished in both eruption series, based on their zoning patterns: type a, a dominant population of unzoned or weakly zoned crystals (Fig. 6a); type b, crystals that are normally zoned from more mag-nesian compositions in the centres, to rims that are the same composition as the unzoned crystals (Fig. 6b); type c, crystals with inner and outer zones of a similar com-position to the unzoned crystals, and an intermediate zone of similar composition to the centres of type b (Fig. 6c). Some crystals of types a and b may be two-dimensional slices of the same three-dimensional crys-tals of type c. We therefore interpret the three crystal types as in general having three zones: margin (zone 1), intermediate (zone 2), and core (zone 3). Unzoned type a

crystals have only zone 1, whereas the zoned type b and c crystals also have zones 2 and 3. The Mg# values of the three zones in the Therasia dacite opx are (from the margin to core) 58–64, 65–71, and 61–65, whereas in the Cape Riva dacite they are 53–58, 60–68, and 53–58.

X-ray element maps of Mg, Fe, Ca and Al of 13 typical opx crystals from the Cape Riva dacite reveal further details of the zoning patterns (Fig. 6). Diffusion rates of these elements decrease in the order Mg–Fe > Ca > Al at about 900°C (Smith & Barron, 1991; Ganguly & Tazzoli, 1994; Klügel, 2001; Nakagawa et al., 2005; Cherniak & Dimanov, 2010; Dohmen et al., 2016), so that whereas Mg–Fe distributions are blurred by subsequent diffusion, Al and Ca distributions provide better images of original growth textures (Tomiya & Takahashi, 2005; Allan et al., 2013). The Al and Ca contents of the different zones were measured by electron microprobe analysis. The weakly zoned crystals commonly exhibit weak oscillatory zonation in Ca and Al (Fig. 6a). One crystal that is normally zoned in Mg–Fe consists of a partially resorbed Al-rich, sector zoned morphology (with oscillatory zoning), infilled and over-grown by Al-poorer opx (Fig. 6d). Comparison of the Al distribution with the maps of Mg and Fe shows that the
<figure><img src="image_8.png" /><figcaption>**Fig. 4.** Histograms of the compositions of plagioclase crystals in the Therasia and Cape Riva dacites. SiO₂ contents in the labels refer to whole-rock compositions presented by Druitt et al. (1999) and Fabbro et al. (2013); n refers to the number of analyses. Crystal cores and rims are grouped together, except for euhedral rim edges in contact with groundmass, which are plotted separately. The Therasia lavas are numbered according to Fabbro et al. (2013). The Cape Riva eruption units (Druitt, 1985) are: A, initial Plinian phase; B, welded ignimbrite; C, nonwelded ignimbrite and lag breccia; D, welded ignimbrite.</figcaption></figure>

sector-zoned core corresponds to the Mg-rich intermediate zone (zone 2) and the Al-poor overgrowth corresponds to the Mg-poorer margin (zone 1); however, subsequent diffusion has smeared out the details of Mg and Fe zoning. Another crystal, which exhibits the three growth zones in Mg-Fe, exhibits the same three zones in Al and Ca

(Fig. 6c). The intermediate zone (zone 2) is richest in Mg, Al and Ca. The outer limit of zone 2 was partially resorbed prior to overgrowth of the margin (zone 1). Despite diffusion of Mg-Fe, a compilation of microprobe analyses shows that the Mg-rich zones are richer in Al (Fig. 7). Ca exhibits a less systematic relationship with Mg and Al.
<figure><img src="image_9.png" /><figcaption>**Fig. 5.** Calibrated BSE images of (a) a plagioclase from the Therasia dacites (GS10-14 XL66) with a calcic core and more sodic rim; (b) a plagioclase from the Cape Riva dacite (GS11-34a XL14) with three sawtooth zones (the outermost one being referred to in the text as the crystal rim). White circles indicate laser analysis spots. The anorthite profile is measured from the calibrated BSE image; trace element profiles were measured by LA-ICP-MS. Continuous lines, measured concentrations; dashed lines, melt equivalent concentrations (see text for explanation of how melt equivalent concentrations were calculated).</figcaption></figure>

# CALCULATION OF MELT COMPOSITIONS

## Plagioclase

### Choice of 'immobile' elements

For the Therasia–Cape Riva series, we are interested in processes that occurred on timescales of less than ~18 kyr (the time between the oldest Therasia dacite and the Cape Riva eruption; Fabbro et al., 2013); elements that remain immobile over this time can be used to

reconstruct melt compositions. We use 'immobile' to mean that the characteristic one-dimensional diffusional length scale, *x*, is less than the size of our 10–15 µm ICP-MS laser-ablation pits, where $x = \sqrt{Dt}$ (*t* is time and *D* is the diffusion coefficient; Fig. 8). The diffusion coefficient strongly depends on temperature, and in the following discussion a temperature of 880°C is used. This is the pre-eruptive temperature of the
<table><caption>**Table 2:** Typical orthopyroxene compositions</caption><thead><tr><th></th><th colspan="3">Therasia (GS10-14 PYX26)</th><th colspan="3">Cape Riva (GS11-34a XL14)</th></tr><tr><th></th><th>Zone 1</th><th>Zone 2</th><th>Zone 3</th><th>Zone 1</th><th>Zone 2</th><th>Zone 3</th></tr></thead><tbody><tr><td colspan="7"><em>Electron microprobe (wt %)</em></td></tr><tr><td>SiO<sub>2</sub></td><td>52.43</td><td>53.41</td><td>52.78</td><td>53.36</td><td>53.38</td><td>52.34</td></tr><tr><td>Al<sub>2</sub>O<sub>3</sub></td><td>1.98</td><td>0.99</td><td>0.66</td><td>0.54</td><td>1.32</td><td>0.52</td></tr><tr><td>FeO</td><td>21.12</td><td>19.55</td><td>23.17</td><td>25.62</td><td>21.99</td><td>24.59</td></tr><tr><td>MgO</td><td>21.72</td><td>23.27</td><td>20.60</td><td>18.11</td><td>21.03</td><td>19.35</td></tr><tr><td>CaO</td><td>1.55</td><td>1.63</td><td>1.84</td><td>1.58</td><td>1.67</td><td>1.59</td></tr><tr><td>TiO<sub>2</sub></td><td>0.46</td><td>0.32</td><td>0.22</td><td>0.20</td><td>0.42</td><td>0.26</td></tr><tr><td>MnO</td><td>0.91</td><td>0.67</td><td>1.04</td><td>1.24</td><td>1.04</td><td>1.23</td></tr><tr><td>Total</td><td>100.49</td><td>99.85</td><td>100.02</td><td>100.63</td><td>100.85</td><td>99.88</td></tr><tr><td>En</td><td>62.2</td><td>65.7</td><td>59.3</td><td>53.9</td><td>60.8</td><td>56.4</td></tr><tr><td>Fs</td><td>34.0</td><td>31.0</td><td>37.4</td><td>42.7</td><td>35.7</td><td>40.2</td></tr><tr><td>Wo</td><td>3.8</td><td>3.3</td><td>3.2</td><td>3.4</td><td>3.5</td><td>3.3</td></tr><tr><td>Mg#</td><td>64.7</td><td>68.0</td><td>61.3</td><td>55.8</td><td>63.0</td><td>58.4</td></tr><tr><td colspan="7"><em>Cations based on 6 oxygens</em></td></tr><tr><td>Si</td><td>1.948</td><td>1.973</td><td>1.979</td><td>2.012</td><td>1.975</td><td>1.986</td></tr><tr><td>Al<sup>iv</sup></td><td>0.052</td><td>0.027</td><td>0.021</td><td>0.000</td><td>0.025</td><td>0.014</td></tr><tr><td>Al<sup>vi</sup></td><td>0.034</td><td>0.016</td><td>0.008</td><td>0.024</td><td>0.033</td><td>0.009</td></tr><tr><td>Fe</td><td>0.656</td><td>0.604</td><td>0.727</td><td>0.808</td><td>0.680</td><td>0.780</td></tr><tr><td>Mg</td><td>1.203</td><td>1.282</td><td>1.152</td><td>1.018</td><td>1.160</td><td>1.094</td></tr><tr><td>Ca</td><td>0.062</td><td>0.065</td><td>0.074</td><td>0.064</td><td>0.066</td><td>0.065</td></tr><tr><td>Ti</td><td>0.013</td><td>0.009</td><td>0.006</td><td>0.006</td><td>0.012</td><td>0.007</td></tr><tr><td>Mn</td><td>0.029</td><td>0.021</td><td>0.033</td><td>0.040</td><td>0.033</td><td>0.040</td></tr></tbody></table>

Therasia and Cape Riva dacites estimated using Fe–Ti oxides and phase equilibria experiments (Fabbro et al., 2013; Cadoux et al., 2014).

Over 18 kyr, La has $x \sim 2$ µm, and can be considered immobile (Cherniak, 2002a). NaSi–CaAl exchange has an $x$ of either 35 nm or 10–15 µm, depending on whether we use the *D* for NaSi–CaAl interdiffusion under dry conditions (Grove et al., 1984) or under water-bearing conditions (Liu & Yund, 1992). In either case, NaSi and CaAl can be considered immobile in plagioclase. The diffusivity of Ti has not been measured, but the high charge of $Ti^{4+}$ ions suggests that it will be slow (Cherniak, 2010). The self-diffusivity of $Si^{4+}$ is similar to that of NaSi–CaAl interdiffusion, and $Ti^{4+}$ may substitute for $Si^{4+}$ in the plagioclase structure (Bindeman et al., 1998; Cherniak, 2003). K and Mg, on the other hand, should have at least partially re-equilibrated with their host melt ($x \sim 700–8000$ µm and 300–5000 µm over 18 kyr, respectively; Giletti & Shanahan, 1997; Van Orman et al., 2014). Strontium exhibits intermediate behaviour; it has $x \sim 50–800$ µm over 18 kyr (Giletti & Casserly, 1994). However, if the near-liquidus residence times for the crystals are significantly less than the $\sim 3$ kyr between the Therasia dome complex and the Cape Riva eruption, then Sr concentrations should be close to growth compositions (Fig. 8a). We conclude that An, La, Ti and possibly Sr are ‘immobile’ on the length and time scales of interest, and (assuming equilibrium element partitioning between crystal and melt) can be used to estimate the composition of the melt from which each plagioclase zone grew. Mg and K, on the other hand, are expected to have undergone partial diffusive re-equilibration.

### Plagioclase–melt trace element partition coefficients

Melt composition cannot be inferred directly from plagioclase An content, which is a function of melt composition, temperature, crystal-growth rate, and $P_{H_2O}$ (Kudo & Weill, 1970; Putirka, 2005). Fluctuations of An may record variations of one or several of these parameters in a manner that is difficult to constrain. On the other hand, melt compositions calculated from plagioclase trace elements (such as Sr, La, Ti, K or Mg) using equilibrium plag–melt partition data are dependent on An content and temperature, but not (it is assumed) on $P_{H_2O}$ (Blundy & Wood, 1991; Bindeman et al., 1998). A global regression of all partition coefficient data by Bédard (2006) suggested that a change of 1–2 wt % $H_2O$ in the melt could produce an $\sim 100$ ppm variation in apparent melt Sr. However, Blundy & Wood (1991) found little difference between the partition coefficients of Sr from experiments both under hydrothermal conditions and with silicate melts, implying that the effect of $P_{H_2O}$ is included in the An term of the partition function. Melt composition may also affect partition coefficients, and this is discussed in more detail for Mg in the Appendix.

The partition coefficients of Bindeman et al. (1998) were used to calculate the Ti, Sr, La, K and Mg contents of melts in equilibrium with measured plagioclase compositions (‘melt-equivalent concentrations’). The partition coefficients ($K$) are described by

$$K = \exp \left( \frac{AX_{An} + B}{RT} \right) \quad (1)$$

where A and B are empirically determined coefficients, $X_{An}$ is the An molar fraction, R is the molar gas constant
<figure><img src="image_10.png" /><figcaption>**Fig. 6.** Mg, Fe, Ca and Al maps of opx phenocrysts from the Cape Riva dacites, with corresponding sketches showing the characteristic features in the Al zoning. White numbers on these sketches are the zone numbers as described in the text. (a) A weakly zoned type *a* crystal (S12-05 PYX26); (b) a type *b* crystal with a partially resorbed high-Al core (S12-05 PYX12); (c) a type *c* crystal with a high-Al mantle (GS11-30b PYX07); (d) a crystal with sector zoning or skeletal growth (S12-06 PYX17). Different types of inclusion can also be distinguished by comparing the different element maps. Glass inclusions (as well as adhering glass) have higher Al, similar Ca and lower Mg and Fe than the opx crystals themselves. Fe–Ti oxides have high Fe, moderately high Al, and low Mg and Ca. Apatite inclusions have high Ca, and low Mg, Fe and Al.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 7.** Compositions of the weakly zoned orthopyroxene crystals, and the different zones of the strongly zoned crystals, from electron microprobe spot analyses. The boundaries between the zones were defined using the Al maps (Fig. 6). All crystals are from the Cape Riva dacite.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 8.** Characteristic diffusion distance against time for species in (a) plagioclase and (b) orthopyroxene, calculated as $x = \sqrt{Dt}$ at 880°C. The upper horizontal dashed lines represent the typical crystal size (1 mm for plag; 500 µm for opx), and the lower horizontal dashed lines represent the size of the analysis (laser pits 10–15 µm in diameter for plag; electron microprobe spot ~1 µm for opx). The left and right vertical dashed lines are the time between the youngest dated Therasia lava and the onset of the Cape Riva eruption (2800 years), and the time between the first dacitic eruption of the Therasia dome complex and the onset of the Cape Riva eruption (18 000 years), respectively. References: 1, NaSi–CaAl interdiffusion at 1500 MPa with 1% H₂O, Liu & Yund (1992); 2, NaSi–CaAl interdiffusion at atmospheric pressure and dry, Grove *et al.* (1984); 3, Cherniak (2002a); 4, Van Orman *et al.* (2014); 5, Giletti & Shanahan (1997); 6, Giletti & Casserly (1994); 7, Ganguly & Tazzoli (1994), along the *a*-axis of orthopyroxene with a Mg# of 60 mol % corrected to the FMQ oxygen buffer using the equation of Allan *et al.* (2013); 8, Dohmen *et al.* (2016), along the *a*-axis of orthopyroxene with a Mg# of 60 mol % at the FMQ oxygen buffer; 9, Nakagawa *et al.* (2005).</figcaption></figure>

and *T* is the temperature in Kelvin. We assumed that *T* was related to An content as follows (Druitt *et al.*, 2012):

$$T = 1128 + 200 \left( \frac{X_{An} - 0.4}{0.4} \right). \quad (2)$$

This implies that plagioclase of An₄₀ is in equilibrium with a melt at 855°C and An₈₀ is in equilibrium at 1055°C, matching the temperatures of Santorini rhyodacites and basaltic andesites, respectively (Michaud *et al.*, 2000). Errors in the calculated melt trace element contents were estimated by standard error propagation, taking into account ICP-MS and electron microprobe counting statistics and published uncertainties on coefficients A and B in the *K* equation of Bindeman *et al.*

(1998). Errors in the magmatic temperatures were estimated by Michaud *et al.* (2000) to be less than about ±20°C, and were not included in the error propagation.

### Interpretation of trace element profiles

The melt-equivalent trace element profiles in Fig. 5 support the conclusion that over the length and time scales of our system Ti, Sr and La are immobile, and Mg is mobile in plagioclase (Fig. 8). The core of the Therasia crystal in Fig. 5a has higher melt Sr, and lower melt La, than the rim. This would be expected from their behaviour in Santorini magmas—there is more Sr in mafic melts and more La in evolved melts (Huijsmans *et al.*, 1988; Druitt *et al.*, 1999). The melt Mg profile, on the other hand, is
<figure><img src="image_13.png" /><figcaption>**Fig. 9.** Melt equivalent trace element contents of plagioclase crystals from the Therasia and Cape Riva eruptions (see text for explanation of how melt equivalent concentrations were calculated). Typical anorthite error bars were calculated from electron microprobe counting statistics. Typical trace element error bars were calculated using LA-ICP-MS counting statistics, estimated errors in the An content, and errors on the partition coefficients given by Bindeman *et al.* (1998). Crystal cores and rims are grouped together, except for euhedral rim edges in contact with groundmass, which are plotted separately.</figcaption></figure>

more complicated. Traced from the margin inwards, it starts correlated with An content; however, the inner part of the core is in equilibrium with a melt with lower Mg than the outer part of the core. This is despite the fact that the inner part of the core has higher An than the outer part. This suggests that Mg in this crystal has partially re-equilibrated. A similar pattern is seen in the crystal in Fig. 5b. The melt Sr values are positively correlated with, and melt La values negatively correlated with, An content. However, in this crystal, the shape of the melt Mg profile also resembles that of An, suggesting that less diffusion has taken place in this crystal compared with Fig. 5a.

## Calculated melt compositions

Melt-equivalent Ti and Sr contents correlate with An (Fig. 9a and b), suggesting that the An content is correlated with the composition and temperature of its host melt, which are themselves also assumed to be correlated through equation (2). Moreover, a plot of melt-equivalent Ti versus melt-equivalent Sr resembles the liquid line of descent of Santorini magmas based on whole-rock compositions (Fig. 10). Whereas Sr

decreases monotonically with increasing differentiation, Ti peaks at intermediate compositions owing to the appearance of Fe–Ti oxides in the liquid line of descent (Druitt *et al.*, 1999). Plagioclase rims in the TH dacites overlap on this plot with whole-rock and groundmass compositions for those same lavas; those of the CR dacite overlap with whole-rock and glass compositions of the Cape Riva pumices, but also extend to more evolved compositions. The data are therefore consistent with derivation of the plagioclase rims from their host melts, with the exception of occasional calcic xenocrysts in the Therasia dacites. These data also suggest that neither Sr nor Ti has diffused significantly since the growth of the crystals.

The An-rich phenocryst cores, however, yield melt Sr contents higher than the whole-rock compositions, indicating derivation from intermediate to mafic melts. The cores are therefore interpreted as either antecrystic or xenocrystic in origin. Those with >280 ppm Sr are probably xenocrystic and were recycled from the Peristeria edifice (530–430 ka), the collapsed remains of which underlie the present-day caldera basin. The Peristeria magmas are the only ones in the last half a
<figure><img src="image_14.png" /><figcaption>**Fig. 10.** Calculated Sr and Ti melt compositions in equilibrium with Therasia and Cape Riva plagioclase crystals. Values calculated from crystal cores and rims are grouped together; those from euhedral rim edges are plotted separately. The shaded areas are the compositional fields of different magma suites from throughout the history of Santorini Volcano. Compositions of Therasia whole-rocks, groundmasses and mafic enclaves, as well as Cape Riva glasses, are taken from Fabbro et al. (2013). The Upper Scoria 2 (US2) data are taken from Druitt et al. (1999), and the Skaros and Peristeria data from Huijsmans (1985). The dashed lines are contours of approximate weight per cent SiO₂. Crystal cores and rims are grouped together, except for euhedral rim edges in contact with groundmass, which are plotted separately.</figcaption></figure>

million years at Santorini to have such high Sr contents (Huijsmans et al., 1988; Druitt et al., 1999). Those with 200–280 ppm Sr lie in the field of the 67–54 ka Skaros shield, onto which the Therasia lavas were discharged, suggesting possible crystal recycling. However, given the similar liquid line of descent in both the Skaros and Therasia lavas, these cores could also be derived from the basaltic parental magmas of the TH dacites.

**Table 3** summarizes the melt compositions from which the cores and rims grew, based on visually cross-correlating An content with melt SiO₂ content via melt Ti content (Fig. 11). The rims of Therasia plagioclase (An₃₉–₅₁) formed in equilibrium with dacitic melts (61–66 wt % SiO₂) of the same composition as the Therasia whole-rock and groundmass. The highest An contents found in cores (An₆₂–₉₁) formed from basaltic to andesitic melts (50–60 wt % SiO₂). Sawtooth zoning within many of these TH crystals records cycling through different melts ranging in composition from andesitic to dacitic. The rims of Cape Riva plagioclase (An₃₆–₄₀) formed in equilibrium with rhyodacitic melts (68–72 wt % SiO₂, equivalent to the host glass), whereas cores (up to An₆₀) record equilibrium with dacitic melts (as low as 61 wt % SiO₂).

## Fingerprinting the origin of plagioclase

Because the major-element compositions of the TH and CR magmas are so similar, we cannot use the major-element compositions of the plagioclase crystals to distinguish between crystals that grew in each magma.

<table><caption>**Table 3:** Estimated melt compositions in equilibrium with the different plagioclase zones</caption><thead><tr><th></th><th>Plagioclase<br>composition<br>(An, mol%)</th><th>Estimated melt<br>composition (SiO<sub>2</sub>, wt%)</th></tr></thead><tbody><tr><td><em>Therasia dacites</em></td><td></td><td></td></tr><tr><td>Calcic cores</td><td>62–91</td><td>50–60 (basalt–andesite)</td></tr><tr><td>Sawtooth peaks</td><td>59–72</td><td>56–64 (andesite–dacite)</td></tr><tr><td>Rims</td><td>39–51</td><td>61–66 (dacite)</td></tr><tr><td><em>Cape Riva dacite</em></td><td></td><td></td></tr><tr><td>Sawtooth peaks</td><td>46–60</td><td>61–67 (dacite)</td></tr><tr><td>Rims</td><td>36–40</td><td>68–72 (rhyodacite)</td></tr></tbody></table>

Sawtooth peaks refer to the most calcic part of the sawtooth zones.

However, differences in incompatible element concentrations of TH and CR dacites do provide a means of differentiating between CR and TH plagioclase. The TH dacites are enriched in incompatible elements (such as La, Ce and K) by 15–25% relative to CR dacite (Fabbro et al., 2013), and this is reflected in the compositions of their plagioclase phenocrysts (Fig. 9c and d). For a given An content, the Therasia plagioclase contains more K and La than the Cape Riva plagioclase, with limited overlap. Sr and Ti contents, on the other hand, show no difference (Fig. 9a and b).

The lower La concentrations in the Cape Riva plagioclase, coupled with the very slow diffusion rates of La (Fig. 8; Cherniak, 2002), indicate that the cores of the Cape Riva plagioclase, like the rims, grew from a CR-like melt, and were not in general inherited from earlier TH magmas. Because K diffuses more rapidly (Fig. 8), some late-stage re-equilibration of this element cannot be ruled out. Only two out of the 90 analyses of 37 Cape Riva plagioclase crystals have sufficiently high K and La to indicate possible derivation from TH magma (Fig. 9c and d), and both of these analyses are from phenocryst cores. The majority of the cores in the CR plagioclase are therefore antecrystic and derived from magma or mushes genetically related to the CR melts.

## Orthopyroxene

The orthopyroxene crystals typically reach 500 µm in length, so that any element with characteristic diffusion length scale $x > 500$ µm will have re-equilibrated at the whole-crystal scale, whereas those with $x < 1$ µm (the size of the electron microprobe beam) will remain close to their initial concentrations. Al has $x \sim 3$ µm over 18 kyr (at 880°C) and can be considered immobile for our purposes (Fig. 8b; Nakagawa et al., 2005). There are no available experimental data for the diffusion coefficient of Ca; however, a value between that of Al and Fe–Mg is assumed by comparison with clinopyroxene data (Cherniak & Dimanov, 2010). Fe and Mg are the fastest diffusing elements in opx considered here, and most opx will have completely re-equilibrated in Fe–Mg over 18 kyr (Allan et al., 2013). However, if the residence times of the crystals are significantly shorter than 18 kyr, then the Fe–Mg contents of the opx may still be
<figure><img src="image_15.png" /><figcaption>**Fig. 11.** (a) Calculated melt Ti concentrations in equilibrium with plagioclase crystals from the Therasia and Cape Riva dacites. (b) Ti concentrations in whole-rock (wr), groundmass separates (gm), and mafic enclaves (en) from the Skaros lavas (Huijsmans, 1985), Upper Scoria 2 eruption (US2; Druitt et al., 1999), Therasia lavas, and Cape Riva eruption (Druitt et al., 1999; Fabbro et al., 2013). (c) and (d) are interpretative sketches of (a) and (b), respectively, with dashed lines showing how plagioclase and whole-rock data can be matched up. Fractionated magmas are those formed primarily through assimilation–fractional crystallization (AFC) processes (Druitt et al., 1999; Fabbro et al., 2013).</figcaption></figure>

close to their original values. For example, $x \sim 160$ µm for Fe–Mg over the ~2800 years between the last Therasia eruption and the Cape Riva eruption, and this drops to $x \sim 30$ µm over 100 years.

Despite the low diffusivity of Al and Ca in opx, it is difficult to use their concentrations to estimate melt compositions in the same way as Sr and Ti in plagioclase. Whole-rock Al contents vary little at Santorini (Huijsmans et al., 1988; Druitt et al., 1999), and the partitioning of Al between the melt and opx is complicated by the fact that Al can occupy both tetrahedral and octahedral sites (Bédard, 2007). Ca is not a trace element in either the crystal or melt, and experimental data for both Al and Ca show that there is no simple relationship between their concentrations in opx and melt composition or other parameters (Bédard, 2007). The Mg# of the opx is therefore more useful as an indicator of melt composition. Figure 12 shows how opx Mg# varies with melt SiO₂ content in phase-equilibria experiments using Cape Riva dacite run by Cadoux et al. (2014). The total range of Mg# will still provide some information about the composition of the melts with which the opx equilibrated. The range of compositions of zones 1, 2 and 3, and the weakly zoned opx, are consistent with equilibration with melts ranging in composition from

the whole-rocks ($\le 67$ wt % SiO₂) to interstitial melts (70–72 wt % SiO₂). This of course assumes equilibrium between opx and melt, which the occurrence of Al sector zoning throws into doubt, as sector zoning involves the precipitation of different opx compositions on different crystal facets owing to local disequilibrium processes (Hollister & Bence, 1967).

## CRYSTAL RESIDENCE TIMESCALES

Residence timescales of plagioclase and opx were calculated at near-liquidus temperatures using the diffusion chronometry method (e.g. Costa et al., 2003, 2008). We distinguish between our calculated 'near-liquidus residence times (NLRT)' and the true residence times of the crystals. The NLRT is a true residence time only if the crystal resided only at the assumed temperature prior to eruption. More complex thermal histories involving prolonged periods at lower temperatures could result in the NLRT underestimating true residence time (Cooper & Kent, 2014). The modelling techniques are described in detail in the Appendix and are summarized below. Typical plagioclase and opx models are shown in Figs 13 and 14, respectively, and the results are summarized in Fig. 15.
<figure><img src="image_16.png" /><figcaption>**Fig. 12.** Mg# of opx crystals plotted against the SiO₂ contents of the coexisting glasses produced in the experiments of Cadoux *et al.* (2014) using the Cape Riva dacite as starting material. The horizontal grey bands represent the compositions of the natural orthopyroxene from the Cape Riva eruption (Druitt *et al.*, 1999; Fabbro *et al.*, 2013). The vertical grey bands are the compositions of the glass from the Cape Riva dacite (Druitt *et al.*, 1999; Fabbro *et al.*, 2013) and the melt compositions calculated from the trace element contents of the most calcic parts of the plagioclase sawtooth zones (Fig. 11, Table 3). The green boxes represent the combinations of glass and opx compositions we would expect if zones 1 and 3 grew from a Cape Riva-like magma whereas zone 2 grew from the magma composition estimated from the sawtooth zones in the plagioclase.</figcaption></figure>

## Mg diffusion in plagioclase

### Modelling

The initial Mg concentrations were estimated using the Ti concentrations of the crystals. Mg and Ti concentrations are correlated in the melt (Druitt et al., 1999; Fabbro et al., 2013), and Ti is in turn correlated with plagioclase An content (Fig. 11). We estimated the expected initial Mg concentration for a given An content using equations (7)–(9) (in the Appendix). The uncertainty in this initial profile was assessed by creating similar equations above and below the best-fit line, such that about two-thirds of the data lay between the two. The equilibrium profiles were calculated using the partition coefficient of Bindeman et al. (1998), assuming that the crystal rims were in equilibrium with the melt [equation (11)]. The Mg diffusion coefficient of Van Orman et al. (2014) was used, as their experiments cover the ranges of temperature and plagioclase An content of Santorini dacites and rhyodacites (850–930°C; Sigurdsson et al., 1990; Fabbro et al., 2013; Cadoux et al., 2014). Where the measured profiles were not completely perpendicular to the crystal rim, we corrected the profile's length using the cosine of the angle between the measured profile and perpendicular. This correction was generally minor (Supplementary Data 1; supplementary data are available for downloading at http://www.petrology.oxfordjournals.org).

The diffusion of Mg was modelled numerically using the method of Costa et al. (2003), implemented using Matlab. The diffusion times were first estimated using one-stage models: the whole crystal was assumed to have grown instantaneously, with no diffusion, then allowed to equilibrate at 880°C. Some of the profiles

could, however, be better fitted by a two-stage model. The two-stage models were run using the technique of Druitt et al. (2012). The core of the crystal was assumed to grow instantaneously, then allowed to partially equilibrate (stage 1) in melt estimated to be in equilibrium with the outermost part of the core. It was then overgrown by the crystal rim with new An and Mg contents, and the whole crystal was allowed to re-equilibrate diffusively with melt of the groundmass composition (stage 2). Diffusion was modelled at temperatures between 880 and 930°C, depending on the An content [equation (2)].

The uncertainty owing to the model fit was assessed by running multiple models for each crystal, and using the shortest and longest times where the modelled profile was still within error of 50% of the measurements. A conservative estimate suggests that the uncertainty in the model fit leads to an uncertainty of a factor of up to 2–4 in the residence time, depending on the crystal. Running models starting with the high and low estimates of the initial Mg concentrations give changes in the NLRT of about ±40 years. Variations of ±25°C lead to changes in NLRT of a factor of two. All these uncertainties were combined to give the ranges of ages shown in Tables 4 and 5.

## Results

Typical one-stage diffusion models are shown in Fig. 13b and e, and a typical two-stage model is shown in Fig. 13c. The results of all the diffusion models, along with calibrated BSE images and An profiles for all the crystals modelled, are presented in Supplementary Data 1. Tables 4 and 5 summarize the results for the Therasia and Cape Riva dacites, respectively.

Most of the measured Mg profiles (especially those in the centres of the crystals) lie between the estimated initial and equilibrium profiles. This confirms that Mg has partially re-equilibrated by diffusion, and allows us to estimate the NLRT. Two crystals from the Cape Riva dacite (Fig. 13e and Supplementary Data Fig. S1.8), however, have Mg concentrations that are within, or close to, the range of uncertainty of the initial Mg concentrations, implying that they resided at high temperature for a very short period prior to eruption.

A common feature of the diffusion models is that one profile appears to be more equilibrated than another in the same crystal (e.g. Supplementary Data Fig. S1.7), and hence gives a longer NLRT. The shorter time is taken as the better estimate, as the models provide the maximum time that the crystals could have resided at a given temperature. Our one-dimensional diffusion models may overestimate the time needed for diffusive re-equilibration (Costa et al., 2003, 2008), as they assume that all diffusion occurs parallel to the modelled profile. Depending on the crystal geometry, and the orientation of the concentration gradients, diffusion in other directions may be significant and may allow the crystal to re-equilibrate faster than the models suggest. Anisotropy in the Mg diffusion coefficient is weak,
<figure><img src="image_17.png" /><figcaption>**Fig. 13.** Diffusion models for two plagioclase crystals from the Cape Riva dacite: (a–c) GS11-30b XL15; (d, e) GS11-34 XL14. (a, d) Anorthite contents along the laser profile, measured using the calibrated BSE image. (b, e) The results of one-stage diffusion models. Dashed lines show how the Mg profile evolves with time. Labels are the time in years for each of the models shown. (c) The results of a two-stage diffusion model. Dashed black line is the modelled result after the first stage (core only); continuous black line is the result after both stages (core plus rim). The vertical dotted line shows the boundary between the core and rim. The core and rim times stated in the diagram are the times for each stage only; the total near-liquidus residence time of the crystals is the core residence time + rim residence time. The model residence times that best constrain the true near-liquidus residence time are highlighted in red.</figcaption></figure>

but may also be enough to explain the difference
(LaTourrette & Wasserburg, 1998; Van Orman et al.,
2014). There may also be other fast paths for diffusion
that are not obvious on BSE images, such as cracks or
grain boundaries. Most notably, the complex geometry

of the core–rim boundary in two-dimensional images
shows that multidimensional effects have probably had
a significant effect on the re-equilibration of the Mg con-
centrations. This means that one-dimensional models
will overestimate the true residence time.
<figure><img src="image_18.png" /><figcaption>**Fig. 14.** Diffusion models for two opx crystals from the Cape Riva dacites. (a–c) S12-06 PYX17; (d, e) GS11-03b PYX08.</figcaption></figure>

Whereas some crystals were best fitted by a one-stage model, most fits were improved by running a two-stage model. For example, in one-stage models of crystal GS11-30b XL15 (Fig. 13b) the calculated profile reaches the measured Mg concentrations at the crystal rim after less than 10 years. However, the Mg concentrations at the crystal centre are not reached until about 450 years. In the two-stage model (Fig. 13c), the crystal core is allowed to equilibrate (without the margin) for 400 years, when it reaches the measured Mg

concentrations (as shown by the dashed line). The crystal rim then grows, and the whole crystal remains in the melt for only a short period (~2 years) before being erupted (as shown by the continuous line). Some Therasia plagioclases (Supplementary Data Figs S1.1, S1.3 and S1.4) can be fitted by two-stage models, where the core first reaches equilibrium with its host melt before being overgrown by the rim.

A few crystals, however, are not well fitted by either one-stage or two-stage models (e.g. GS10-14 XL32,
<figure><img src="image_19.png" /><figcaption>**Fig. 15.** Near-liquidus residence timescales of plagioclase and opx crystals from (a, b) the Therasia and (c, d) Cape Riva dacites. (a, c) Opx timescales; (b, d) plagioclase timescales. Plagioclase rim residence times are either the duration of the second stage of the two-stage models or the residence times calculated from one-stage models where two-stage models were not run. Plagioclase core residence times are the combined durations of the first and second stages of the two-stage models. Error bars represent the uncertainties on the NLRT based on uncertainties in the model fit, the initial conditions, and the temperature (as discussed in the text). The average time between Therasia eruptions, and the time between the last Therasia eruption and the Cape Riva eruption, was calculated using radiometric dates (Fabbro et al., 2013). Opx residence timescales are calculated using the diffusion coefficients of Allan et al. (2013, shorter NLRT) and Dohmen et al. (2016, longer NLRT).</figcaption></figure>

Supplementary Data Fig. S1.3), which again may be because not all the diffusion occurred parallel to the measured profile. For GS10-14 XL32, it is possible that there is a boundary between two crystals at about 250 µm from the crystal edge along profile 1, and this would provide a fast path for Mg. Another possibility is that the crystal grew in multiple stages rather than just two. Both explanations would require us to reduce the estimate of the NLRT.

Repeated cycles of crystal growth, partial dissolution, re-equilibration, and renewed growth will affect the shape of the Mg profile, and the one- and two-stage models presented here are clearly simplifications. However, modelling the diffusion of Mg in plagioclase does provide important constraints on the timescales of magmatic processes recorded by the zoning patterns of

the crystals. Most of these processes will cause our models to overestimate NLRT. The NLRT presented here, with the exception of the few cores that equilibrated totally with their host melts, must be considered maximum estimates of the true residence time.

### Calculated plagioclase timescales

The resulting NLRT values for the Therasia and Cape Riva plagioclase are compiled in Tables 4 and 5, and shown in Fig. 15. We focus here on the timescales from the two-stage models, distinguishing the NLRT given by the crystal core (total time of two-step models) from that given by the crystal rim (second stage of two-step models).

Let us reflect on the geological meaning of these two timescales. Because the antecrystic or xenocrystic cores
<table><caption>**Table 4:** Near-liquidus residence times, in years, calculated in diffusion models of Mg in plagioclase crystals from the Therasia dacites</caption><thead><tr><th rowspan="2">Crystal</th><th rowspan="2">Profile</th><th colspan="2">1 Stage</th><th colspan="2">2 Stage</th><th colspan="2"></th></tr><tr><th></th><th></th><th>Core</th><th></th><th>Rim</th><th></th></tr></thead><tbody><tr><td rowspan="3">Flow 5 (GS10-17)</td><td>XL60</td><td>—</td><td>2200</td><td>(630-9000)</td><td><strong>3000</strong></td><td>(>1000)</td><td><strong>15</strong></td><td>(0-120)</td></tr><tr><td>1</td><td></td><td>800</td><td>(100-3400)</td><td></td><td></td><td></td><td></td></tr><tr><td>2</td><td></td><td>330</td><td>(10-1100)</td><td><strong>125</strong></td><td>(30-480)</td><td><strong>50</strong></td><td>(0-330)</td></tr><tr><td rowspan="3">Flow 8 (GS10-14)</td><td>XL32</td><td>1</td><td>1100</td><td>(200-3700)</td><td>500</td><td>(>150)</td><td>75</td><td>(5-330)</td></tr><tr><td>2</td><td></td><td>800</td><td>(25-2100)</td><td><strong>500</strong></td><td>(>7)</td><td><strong>75</strong></td><td>(5-480)</td></tr><tr><td>XL66</td><td>1</td><td>9000</td><td>(>430)</td><td>5000</td><td>(>1200)</td><td>400</td><td>(130-1100)</td></tr><tr><td></td><td></td><td>2</td><td>3000</td><td>(>25)</td><td><strong>1500</strong></td><td>(>400)</td><td><strong>15</strong></td><td>(0-130)</td></tr><tr><td>Flow 9 (GS10-16)</td><td>XL148</td><td>—</td><td><strong>180</strong></td><td>(10-1700)</td><td></td><td></td><td></td><td></td></tr></tbody></table>

Models in bold are those considered to best estimate the NLRT of each crystal. The range of possible ages taking into account the uncertainties in the text is shown in parentheses.

<table><caption>**Table 5:** Near-liquidus residence times, in years, calculated in diffusion models of Mg in plagioclase crystals from the Cape Riva dacite</caption><thead><tr><th rowspan="2">Crystal</th><th rowspan="2">Profile</th><th colspan="2">1 Stage</th><th colspan="2">2 Stage</th><th colspan="2"></th></tr><tr><th></th><th></th><th>Core</th><th></th><th>Rim</th><th></th></tr></thead><tbody><tr><td rowspan="2">CR-A (GS11-34a)</td><td>XL09</td><td>1</td><td>20</td><td>(5-240)</td><td><strong>10</strong></td><td>(2-120)</td><td><strong>5</strong></td><td>(0-100)</td></tr><tr><td></td><td>2</td><td>70</td><td>(10-400)</td><td>15</td><td>(5-120)</td><td>10</td><td>(0-110)</td></tr><tr><td rowspan="2">XL14</td><td></td><td>1</td><td><strong>0</strong></td><td>(0-120)</td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>2</td><td>80</td><td>(15-500)</td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="2">XL23</td><td></td><td>1</td><td>45</td><td>(0-260)</td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>2</td><td><strong>35</strong></td><td>(0-400)</td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="2">CR-B (S12-06)</td><td>XL24</td><td>1</td><td>300</td><td>(45-1500)</td><td>10</td><td>(>5)</td><td>250</td><td>(55-680)</td></tr><tr><td></td><td>2</td><td>170</td><td>(20-480)</td><td><strong>125</strong></td><td>(30-380)</td><td><strong>1</strong></td><td>(0-4)</td></tr><tr><td rowspan="2">CR-C (GS11-30b)</td><td>XL15</td><td>—</td><td>180</td><td>(5-980)</td><td><strong>400</strong></td><td>(130-1100)</td><td><strong>2</strong></td><td>(1-16)</td></tr><tr><td>XL18</td><td>—</td><td>20</td><td>(5-280)</td><td><strong>50</strong></td><td>(5-280)</td><td><strong>25</strong></td><td>(5-140)</td></tr><tr><td rowspan="2">CR-C (S12-05)</td><td>XL03</td><td>1</td><td>25</td><td>(2-180)</td><td><strong>10</strong></td><td>(4-180)</td><td><strong>4</strong></td><td>(0-12)</td></tr><tr><td></td><td>2</td><td>220</td><td>(20-1100)</td><td>75</td><td>(5-280)</td><td>1</td><td>(0-16)</td></tr><tr><td rowspan="2">CR-D (GS11-39b)</td><td>XL29</td><td>1</td><td>60</td><td>(10-680)</td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>2</td><td><strong>4</strong></td><td>(1-22)</td><td></td><td></td><td></td><td></td></tr></tbody></table>

Models in bold are those considered the best estimate the NLRT of each crystal. The range of possible ages taking into account the uncertainties in the text is given in parentheses.

could have had protracted histories at variable temperatures, we can attach no meaning to the derived stage 1 timescale in general other than that it provides an upper limit to the total length of time that the core remained at near-liquidus temperatures (prior to, and following, rim overgrowth) before eruption. The significance of the stage 1 timescales may also vary from crystal to crystal. Some of the cores are clearly recycled, and the timescales recovered from these may record multiple episodes of storage and remobilisation. However, some of the cores originated within similar magmas to their hosts, and the stage 1 timescales from these crystals may reflect events associated with the final assembly of the erupted melt body. However, the stage 2 timescales record the growth and subsequent near-liquidus residence of the rims, and provide the clearest record of pre-eruptive processes.

In the Therasia plagioclase, the cores yield NLRT estimates of 125-3000 years (total range with uncertainties is 30 to >3000 years), whereas the rims give only 15-75 years (total range 0-480 years). The one-stage model timescales are consistent with the core timescales from two-stage modelling. In the Cape Riva plagioclase, timescales given by the crystal cores are 10-400 years (2-1100 years), and those from the rims are 1-25 years (0-140 years).

## Mg-Fe diffusion in orthopyroxene

### Modelling

Fe-Mg concentrations were modelled in opx across individual zone boundaries. The initial zone boundaries were assumed to be step functions. The initial Mg# on either side of the boundary was estimated by taking the values of Mg# at a distance away from the boundary, where they appeared to be unmodified by diffusion.
There are two formulations for the diffusion coefficient of Mg in opx: that of Allan et al. (2013), which is based on the observations of Ganguly & Tazzoli (1994) and Klügel (2001), and one based upon more recent experiments by Dohmen et al. (2016). If we use the formulation of Allan et al. (2013), we need to extrapolate from very different oxygen fugacities, as the original experiments of Ganguly & Tazzoli (1994) were carried out at the iron-wüstite (IW) buffer. However, Dohmen et al. (2016) used only Mg-rich opx (Mg# = 91–98). As both the oxygen fugacity and the Fe content of the opx have a large impact on the diffusivity (Cherniak & Dimanov, 2010), it is unclear to us which formulation is more accurate for our crystals. Most of the NLRT discussed here use the diffusion coefficient of Allan et al. (2013); however, we also present results using the diffusion coefficient of Dohmen et al. (2016) for comparison. NLRT were calculated at 880°C and at the fayalite–magnetite–quartz (FMQ) oxygen fugacity buffer (Fabbro et al., 2013; Cadoux et al., 2014). Because the variation in Mg# in each profile was always small (<10%), the dependence of the Mg-Fe diffusion coefficient on Mg# within a single profile was ignored, and we used a value halfway between the Mg# on either side of the step function. This means that the profiles can be fitted by a simple analytical solution to the diffusion equation (Crank, 1957; Lasaga, 1998). Profiles perpendicular to zone boundaries were measured on calibrated BSE images using the ImageJ software (Rasband, 1997), and were fitted by eye in a spreadsheet program.

The uncertainty owing to the model fit was assessed in a similar fashion to the plagioclase crystals: using the largest range in times over which the modelled profile was still within the range of scatter in the data for at least 50% of the boundary. The uncertainty in the NLRT owing to temperature and oxygen fugacity was calculated directly from their effects on the diffusion coefficient. Using a temperature of 855°C with the diffusion coefficient of Allan et al. (2013) increases the NLRT by a factor of 1.4. Increasing the temperature to 905°C reduces the NLRT by a factor of 0.7. Changes in oxygen fugacity of ±0.3 log units change the NLRT by ±12%. The diffusion coefficient calculated by Dohmen et al. (2016) is more sensitive to temperature, but less sensitive to oxygen fugacity: changing the temperature by ±25°C changes the NLRT by a factor of two, whereas changing the oxygen fugacity by ±0.3 log units changes the NLRT by ±4%.

## Results

Most of the Mg-Fe profiles can be modelled as initial step functions that have broadened by diffusion (Fig. 14 and Supplementary Data 2). However, there is evidence that not all the gradients are entirely due to diffusion. Our approach in constructing models for these profiles was to estimate the maximum possible time that the observed compositional gradients could have survived at near-liquidus temperatures (880°C). Assuming that

the gradient was produced entirely by diffusion in one dimension through an infinitely long solid, starting from an initial step function, provides a robust maximum estimate of the NLRT, even if one or more of those assumptions turns out to be incorrect.

There are occasionally large differences in the times calculated for profiles along different directions in the same crystal (e.g. Fig. 14a–c), which could be due to anisotropy in diffusion or growth rates. Fe-Mg diffusion in opx is faster along the *c*-axis than the *a*- and *b*-axes (Ganguly & Tazzoli, 1994). Although the magnitude of this anisotropy has not been well constrained, Schwandt et al. (1998) found that the self-diffusion coefficient of Mg varied by less than 0.5 log units below 900°C, and Dohmen et al. (2016) found that it varied by a factor of 3.5 (0.54 log units) between 870 and 1100°C. Diffusion anisotropy could explain some of the variations in apparent NLRT seen within single crystals, but not all. However, the models assume that the Mg# profile was initially a step function, so that any growth zoning will increase the apparent NLRT. For example, longer NLRT given by models along the long axes of the opx can be explained by faster growth in that direction. Comparison of the Fe-Mg zoning patterns with those of Al support this possibility. Al diffuses so slowly in opx that its concentrations will reflect growth conditions only, and the Al zone boundaries are generally more diffuse in profiles along the long dimension of the crystal (e.g. Fig. 6d). Using the shape of the Al profiles to estimate the initial Fe-Mg profile, we find that growth zoning could completely explain the time differences between profiles along different directions in the same crystal. However, the uncertainties involved—particularly owing to the low resolution of the Al maps with respect to the size of the diffusion profiles—mean that some diffusion anisotropy cannot be ruled out. This is consistent with the results of Allan et al. (2013), who also found that gradients along the *c*-axis (parallel to the direction of elongation of the opx crystals) were best explained by a mix of growth and diffusion. Where only one profile in a crystal was modelled, it was chosen to be perpendicular to the long axis to minimize this effect. Choosing such a profile, far away from crystal corners, also reduced the impact of diffusion in three dimensions, perpendicular to the profile (Krimer & Costa, 2017). Where more than one profile was modelled, the shorter residence time was used as the best estimate of NLRT. The sharpness of the zoning patterns in most of the Al maps of the strongly zoned crystals suggests that the effect of growth on the timescales calculated for these crystals is small. This is less true for the weakly zoned crystals, as the Al and Ca often show gradual changes in concentration. However, correcting for the effect of growth will always reduce the NLRT, so the values presented here are maxima.

The NLRT of the weakly zoned type *a* crystals are more difficult to interpret than those of the type *b* and *c* crystals. The zoning patterns of the type *a* crystals are irregular, and do not appear to represent changes in the
local environment during growth. The opx do have many apatite and magnetite inclusions that are usually surrounded by small pockets of trapped melt, and the weak zoning in the type *a* crystals may be due to the re-equilibration of these inclusions with the host opx. It is also possible that the type *a* crystals were originally more strongly zoned, and they have partially re-equilibrated, erasing any evidence of their previous compositions. However, there are two lines of evidence that suggest that they were never strongly zoned. First, all the Mg-rich zones of the strongly zoned type *b* and *c* crystals also have high Al. The diffusivity of Al is too low for it to have been significantly altered (Nakagawa *et al.*, 2005), yet we do not see Al-rich zones in the type *a* crystals. Second, the NLRT from the type *a* crystals matches those from the type *b* crystals, which suggests that the type *a* and *b* crystals are recording the same events.

There are a number of features in the measured profiles that demonstrate that our models do not capture all the complexity of the Fe–Mg zoning. For example, the NLRT estimated for the inner, zone 3–zone 2 boundary of many crystals are shorter than those estimated for their outer, zone 2–zone 1 boundary. The boundary between zone 2 and zone 1 also occurs in two steps in some crystals. Finally, diffusion was modelled at a constant temperature of 880°C, whereas above it was argued that the Mg-rich zone 2 grew from a less evolved melt than zones 1 and 3. The implications of each of these complications are discussed in more detail in the Appendix; however, in each case the effect is the same: it leads to our models overestimating the NLRT.

### Calculated opx timescales

The opx NLRT are presented in Tables 6 and 7, and plotted alongside plagioclase NLRT in Fig. 15. For both Therasia and Cape Riva opx, we plot NLRT given by (1) the boundary between zones 3 and 2 in zoned crystals (12 crystals), (2) the boundary between zones 2 and 1 in zoned crystals (28 crystals), and (3) some minor zone boundaries within weakly zoned phenocrysts (13 crystals). Using the diffusion coefficient of Allan *et al.* (2013), values for both boundaries in zoned Therasia opx range from <1 to 30 years (total range with uncertainties <1–90 years), whereas those for weakly zoned opx and both boundaries in zoned opx in Cape Riva pumices range from <1 to 40 years (<1–120 years).

In the Cape Riva dacite, the type *c* crystals (those with all three zones) have shorter NLRT than the type *a* and *b* crystals. Both the inner and outer boundaries of the type *c* crystals have NLRT of 0.09–5 years (0.03–11 years). There are occasional type *a* and *b* crystals that have NLRT similar to the type *c* crystals—these may be slices through the outer zones of a type *c* crystals—however, most have NLRT of 5–40 years (1.2–220 years).

The choice of diffusion coefficient makes a large difference to the results, with the diffusion coefficient of Dohmen *et al.* (2016) giving NLRT 30–40 times greater

<table><caption>**Table 6:** Near-liquidus residence times, in years, calculated in diffusion models of Fe–Mg in orthopyroxene from the Therasia dacites</caption><thead><tr><th>Crystal</th><th>Profile</th><th>Boundary between zones 3 and 2</th><th>Boundary between zones 2 and 1</th></tr></thead><tbody><tr><td colspan="4"><strong>Flow 5 (GS10-17)</strong></td></tr><tr><td>PYX22</td><td>—</td><td>2.6 (0.9-7)</td><td>2.4 (0.5-7)</td></tr><tr><td>PYX24</td><td>—</td><td>—</td><td>24 (4-60)</td></tr><tr><td>PYX30</td><td>1</td><td>—</td><td>6 (1.8-60)</td></tr><tr><td></td><td>2</td><td>—</td><td>14 (5-40)</td></tr><tr><td colspan="4"><strong>Flow 8 (GS10-14)</strong></td></tr><tr><td>PYX01</td><td>—</td><td>—</td><td>6 (1.2-14)</td></tr><tr><td>PYX19</td><td>—</td><td>—</td><td>1.1 (0.3-2.8)</td></tr><tr><td>PYX24</td><td>—</td><td>—</td><td>2.5 (0.6-8)</td></tr><tr><td>PYX26</td><td>—</td><td>0.17 (0.08-0.4)</td><td>0.27 (0.06-0.9)</td></tr><tr><td>PYX29</td><td>—</td><td>—</td><td>20 (8-60)</td></tr><tr><td colspan="4"><strong>Flow 9 (GS10-16)</strong></td></tr><tr><td>PYX02</td><td>1</td><td>—</td><td>30 (11-90)</td></tr><tr><td></td><td>2</td><td>—</td><td>90 (30-190)</td></tr><tr><td>PYX03</td><td>1</td><td>—</td><td>7 (2.2-20)</td></tr><tr><td></td><td>2</td><td>—</td><td>10 (3-27)</td></tr><tr><td>PYX05</td><td>—</td><td>3 (1.8-7)</td><td>6 (1.8-16)</td></tr></tbody></table>

In crystals with multiple profiles, only the profile that provides the best estimate of the residence time (as discussed in the text) is used. The other profiles in those crystals are presented here in italics, but are not considered further.

than the diffusion coefficient of Allan *et al.* (2013). Using Dohmen *et al.* (2016), the range of NLRT of the Therasia opx crystals becomes 6–1000 years (2–4000 years), and the range of NLRT of the Cape Riva opx crystals becomes 4–1700 years (1–13 000 years).

### Comparison of plagioclase and orthopyroxene timescales

Near-liquidus residence timescales for the plagioclase rims and opx [diffusion coefficient of Allan *et al.* (2013)] in both Therasia and Cape Riva dacites range from less than a year to a few decades, possibly up to a couple of centuries if the maximum error is assumed. The plagioclase rim NLRT are similar to those of type *a* and *b* opx although the type *c* opx have shorter NLRT. Shorter times for opx could be explained as an artefact of the tabular shape and complex resorption surfaces found in the plagioclase, which will make diffusion in two and three dimensions significant and increase the true rate of diffusion relative to our 1D models (e.g. Supplementary Data Fig. S1.4; Costa *et al.*, 2003, 2008). In contrast, the opx crystals have more prismatic shapes and planar zone boundaries, which make the use of one-dimensional diffusion more realistic. The higher resolution of our opx transects also allows us to more precisely constrain shorter timescales; the 10–15 µm size of the laser pits in the plagioclase crystals makes it difficult to resolve NLRT of less than a year.

In contrast, opx NLRT using the diffusion coefficient of Dohmen *et al.* (2016) are significantly longer than the plagioclase rim NLRT. This could be explained in three different ways: (1) the opx grew before the plagioclase rims; (2) the opx models overestimate the true NLRT by more than the plagioclase models do; or (3) at least one
<table><caption>**Table 7:** Near-liquidus residence times, in years, calculated in diffusion models of Fe–Mg in orthopyroxene from the Cape Riva dacite</caption><thead><tr><th>Crystal</th><th>Profile</th><th>Boundary between zones 3 and 2</th><th>Boundary between zones 2 and 1</th><th>Weakly zoned</th></tr></thead><tbody><tr><td colspan="5"><strong>Phase A (GS11-34a)</strong></td></tr><tr><td>PYX36</td><td>—</td><td>—</td><td>—</td><td>1.2 (0.28–3)</td></tr><tr><td>PYX39</td><td>1</td><td>—</td><td>9 (3–21)</td><td>—</td></tr><tr><td>PYX39</td><td>2</td><td>—</td><td>250 (100–500)</td><td>—</td></tr><tr><td>PYX59</td><td>1</td><td>1.5 (0.4–4)</td><td>0.5 (0.12–1.5)</td><td>—</td></tr><tr><td>PYX59</td><td>2</td><td>15 (2.7–30)</td><td>1.5 (0.3–8)</td><td>—</td></tr><tr><td>PYX62</td><td>—</td><td>—</td><td>40 (9–220)</td><td>—</td></tr><tr><td>PYX63</td><td>1</td><td>5 (1.5–11)</td><td>4 (1.2–10)</td><td>—</td></tr><tr><td>PYX63</td><td>2</td><td>0.6 (0.21–1.5)</td><td>0.09 (0.03–0.5)</td><td>—</td></tr><tr><td>PYX64</td><td>—</td><td>—</td><td>—</td><td>40 (13–120)</td></tr><tr><td>PYX65</td><td>—</td><td>—</td><td>—</td><td>30 (10–70)</td></tr><tr><td colspan="5"><strong>Phase B (S12-06)</strong></td></tr><tr><td>PYX01</td><td>1</td><td>—</td><td>5 (1.2–11)</td><td>—</td></tr><tr><td>PYX01</td><td>2</td><td>—</td><td>1.7 (0.6–6)</td><td>—</td></tr><tr><td>PYX04</td><td>—</td><td>—</td><td>5 (1.2–10)</td><td>—</td></tr><tr><td>PYX07</td><td>—</td><td>—</td><td>20 (6–50)</td><td>—</td></tr><tr><td>PYX13</td><td>—</td><td>2.7 (1.0–7)</td><td>—</td><td>—</td></tr><tr><td>PYX17</td><td>1</td><td>—</td><td>11 (4–24)</td><td>—</td></tr><tr><td>PYX17</td><td>2</td><td>—</td><td>30 (10–60)</td><td>—</td></tr><tr><td>PYX17</td><td>3</td><td>—</td><td>450 (160–940)</td><td>—</td></tr><tr><td colspan="5"><strong>Phase C (GS11-30b)</strong></td></tr><tr><td>PYX01</td><td>—</td><td>—</td><td>13 (4–28)</td><td>—</td></tr><tr><td>PYX03</td><td>—</td><td>—</td><td>—</td><td>8 (1.4–30)</td></tr><tr><td>PYX07</td><td>—</td><td>—</td><td>—</td><td>14 (2.8–40)</td></tr><tr><td>PYX08</td><td>1</td><td>6 (2.4–14)</td><td>8 (2.0–30)</td><td>—</td></tr><tr><td>PYX08</td><td>2</td><td>0.22 (0.08–0.5)</td><td>0.23 (0.04–0.6)</td><td>—</td></tr><tr><td>PYX10</td><td>1</td><td>—</td><td>—</td><td>6 (1.2–16)</td></tr><tr><td>PYX10</td><td>2</td><td>—</td><td>—</td><td>9 (2.6–22)</td></tr><tr><td>PYX29</td><td>—</td><td>—</td><td>30 (10–60)</td><td>—</td></tr><tr><td colspan="5"><strong>Phase C (S12-05)</strong></td></tr><tr><td>PYX01</td><td>1</td><td>0.7 (0.20–1.8)</td><td>1.3 (0.22–3)</td><td>—</td></tr><tr><td>PYX01</td><td>2</td><td>1.8 (0.6–3)</td><td>5 (1.8–16)</td><td>—</td></tr><tr><td>PYX06</td><td>1</td><td>—</td><td>—</td><td>2.9 (0.8–7)</td></tr><tr><td>PYX06</td><td>2</td><td>—</td><td>—</td><td>14 (4–30)</td></tr><tr><td>PYX12</td><td>—</td><td>—</td><td>22 (5–50)</td><td>—</td></tr><tr><td>PYX18</td><td>1</td><td>1.7 (0.22–3)</td><td>2.1 (0.8–6)</td><td>—</td></tr><tr><td>PYX18</td><td>2</td><td>0.16 (0.05–0.30)</td><td>0.22 (0.05–0.7)</td><td>—</td></tr><tr><td>PYX17</td><td>—</td><td>—</td><td>—</td><td>5 (2–20)</td></tr><tr><td>PYX26</td><td>1</td><td>—</td><td>—</td><td>19 (6–60)</td></tr><tr><td>PYX26</td><td>2</td><td>—</td><td>—</td><td>11 (4–40)</td></tr><tr><td>PYX28</td><td>1</td><td>0.9 (0.22–2.1)</td><td>4 (1.0–13)</td><td>—</td></tr><tr><td>PYX28</td><td>2</td><td>—</td><td>0.9 (0.16–3)</td><td>—</td></tr><tr><td>PYX34</td><td>—</td><td>—</td><td>—</td><td>4 (0.6–10)</td></tr><tr><td colspan="5"><strong>Phase D (GS11-39b)</strong></td></tr><tr><td>PYX11</td><td>—</td><td>—</td><td>—</td><td>19 (6–40)</td></tr><tr><td>PYX12</td><td>—</td><td>—</td><td>—</td><td>11 (2.8–30)</td></tr><tr><td>PYX21</td><td>—</td><td>5 (1.2–11)</td><td>—</td><td>—</td></tr><tr><td>PYX31</td><td>—</td><td>0.16 (0.06–0.5)</td><td>1.5 (0.5–2.8)</td><td>—</td></tr><tr><td>PYX37</td><td>—</td><td>—</td><td>12 (5–24)</td><td>—</td></tr></tbody></table>

In crystals with multiple profiles, only the profile that provides the best estimate of the residence time (as discussed in the text) is used. The other profiles in those crystals are presented here in italics, but are not considered further.

of the formulations for the diffusion coefficients we used does not apply to our system. There is no textural evidence of the opx crystallizing before the plagioclase rims, and given the uncertainties in the diffusion coefficients of both Mg in plagioclase and Fe–Mg in opx and in our estimates of the NLRT, we infer that the plagioclase rims and opx crystallized concomitantly over similar timescales.

Plagioclase antecrystic or xenocrystic cores give NLRT of centuries to millennia. However, if the cores had spent even some of that time at lower temperatures

(in a more crystalline environment), then the true residence times would be much longer.

## DISCUSSION

### Processes recorded by crystal textures and compositions

The cores of the plagioclase crystals have had long and complex histories, as shown by their textural (multiple sawtooth zones, internal resorption surfaces, patchy
zoning, An-rich inner cores), and compositional diver-
sities. They also give model residence timescales of up
to ~10³ years, even when calculated at near-liquidus
temperatures, suggesting long pre-eruptive storage.
Most of the cores are antecrystic, in that their K and La
contents clearly link them genetically to their host
magma series; however, a few cores in TH dacites have
high Sr contents indicative of derivation from the old
Peristeria (530-430 ka) volcano, and are probably
xenocrystic.

Our study offers few firm constraints on the origins
of the crystal cores. The variation of An content appears
to be mostly controlled by melt composition and tem-
perature, rather than P<sub>H₂O</sub> [as also found by Cashman &
Blundy (2013) for Mt St Helens], because there is a
strong correlation between An and calculated melt Sr
(Fig. 9a) and the latter is in principle little affected by
P<sub>H₂O</sub>. Cores in the Therasia plagioclase grew from mafic
andesite to dacite melts, and those in the Cape Riva
from silicic andesite to dacite. There is little evidence for
storage of the cores in low-temperature, high-crystallin-
ity environments. For example, there are no inclusions
of low-temperature accessory minerals such as zircon,
and there is no zoning to more evolved, incompatible
element-enriched compositions near the edge of the
cores. However, some cores may have been derived
from gabbroic-dioritic mush bodies with perhaps 30-
60% melt from which the TH and CR silicic melts were
probably extracted (Mahood, 1990; Sparks et al., 1990;
Koyaguchi & Kaneko, 1999, 2000; Hildreth & Wilson,
2007; Huber et al., 2009; Dufek & Bachmann, 2010;
Cooper & Kent, 2014). If so, the mush zones were com-
positionally distinct, because the La and K contents
show that cores in the CR plagioclases crystallized from
CR-type (incompatible-poorer) magmas, not TH-type
(incompatible-richer) magmas. Calculated NLRTs for
the cores could reflect the duration of residence inte-
grated over the entire temperature-variable history of
the mush, so that the actual residence times could be
much longer than those calculated at near-liquidus tem-
peratures in Tables 4 and 5 (Ruprecht & Cooper, 2012;
Cooper & Kent, 2014). Another possibility is that the
cores grew during prolonged cycling through melt-
dominated volumes of different compositions and
temperatures.

The plagioclase rims record the last phase of crystal-
lization prior to eruption. The normal zoning of the rims
is consistent with crystallization from the host magmas,
suggesting that the opx and plagioclase rims grew sim-
ultaneously. This is supported by the similarities in
NLRT for the opx [<1-30 years for Therasia and <1-
40 years for Cape Riva, using Allan et al. (2013)] and for
the plagioclase rims (15-180 years for Therasia; <1-
35 years for Cape Riva). It seems more likely that the
compositional variations within the phenocryst rims
may record the cycling of these crystals through melt-
dominated magma reservoirs of different compositions
during growth, perhaps as different magma batches
were assembled into the final pre-eruptive reservoir.

Zoned plagioclase and opx crystals reveal the involve-
ment of diverse melt compositions (andesite to dacite
for the Therasia magmas; dacite to rhyodacite for the
Cape Riva magma). However, the large number of unz-
oned crystals suggests that the volume of the more
mafic magma was smaller than that of the more
evolved magma, as these crystals 'saw' only the more
evolved compositions during their growth. The NLRT
for opx using the diffusion coefficient of Dohmen et al.
(2016) are longer (4-1700 years), but are still within
error of the NLRT of the plagioclase rims. Furthermore,
opx in both dacite types appears to have grown rapidly.
Sector-zoned morphologies are visible in Al compos-
itional images from the CR dacite (Fig. 6d), and labora-
tory experiments have shown that the production of Al
sector zoning in opx occurs only during rapid growth
(Schwandt & McKay, 2006). The NLRT using the diffu-
sion coefficient of Allan et al. (2013) are consistent with
rapid opx growth followed by monthly, yearly or dec-
adal residence at high temperature.

Several lines of evidence suggest that the Therasia
dacites were fed by multiple small magma reservoirs,
rather than a single, long-lived reservoir. First, whole-
rock and melt compositions vary slightly throughout
the Therasia sequence, along with temperature and
oxygen fugacity, but not in any systematic way (Fabbro
et al., 2013). Second, the compositions of the euhedral
plagioclase rims vary from one dacite flow to another
(Fig. 4). Third, the plagioclase crystals have distinctive
textures in different lava flows. For example, in some
flows many crystals contain large, brown glass inclu-
sions, but in other flows these inclusion-rich plagioclase
crystals are absent. Plagioclase crystals with sawtooth-
zoned cores are also present in only some of the lavas.
We therefore see little crystal recycling between the dif-
ferent Therasia lavas that could indicate that they
shared a common storage region. The Therasia erup-
tions appear instead to have been fed by the ascent of
multiple volumes of dacite, each containing a slightly
different cargo of crystals.

In contrast, the Cape Riva eruption took place from a
single well-mixed magma reservoir. The whole-rock
and glass compositions of samples taken from through-
out the eruptive sequence are uniform to within analyt-
ical error, excluding the hybrid andesite that occurs as a
minor component in the Plinian phase products and
possibly records a late-stage recharge event (Druitt
et al., 1999; Fabbro et al., 2013). The plagioclase pheno-
crysts are all texturally similar and have composition-
ally restricted rim compositions that are independent of
eruptive phase.

The processes that drove phenocryst growth

A fundamental observation of our study is that crystals
from the Therasia eruptions (multiple discharges of
small volumes of lava over 15kyr) and those from the
much larger, caldera-forming Cape Riva eruption show
many similarities: (1) the plagioclases are similar
texturally and compositionally; (2) the opx are also simi-
lar texturally and compositionally; (3) plagioclase rims
and opx from both systems give mostly annual to dec-
adal timescales of phenocryst residence. This implies
intrinsic processes of phenocryst growth that were
common to both systems. What mechanism of crystal
growth common to these two events of very different
magnitude can explain these similarities?

Crystal growth in magmas can be driven by a de-
crease in temperature or by decompression under H₂O-
saturated conditions. Santorini melts are known to be
saturated with an H₂O-rich vapour phase in the upper
crust (Cadoux et al. 2014; Druitt et al., 2016), so that
vapour-saturated magma ascent and outgassing is to
be expected.

Recent work combining phase equilibria studies of Santorini magmas (a basalt, an andesite, a dacite and three rhyodacites; Cadoux et al., 2014; Andújar et al., 2015, 2016) and H₂O + CO₂ barometry of crystal-hosted melt inclusions (Cadoux et al., 2014; Druitt et al., 2016) has enabled us to place constraints on the generation and pre-eruptive storage of silicic magmas beneath the volcano. Mantle-derived basalt intrudes at about 400 MPa (16 km using a density of 2600 kg m⁻³), at the interface between the lower and upper crust, where it fractionates to intermediate compositions. Silicic melts are generated in the 400–250 MPa pressure range (16–10 km), then transferred to shallow pre-eruptive storage reservoirs (4–6 km) prior to Plinian eruptions.

We can use the phase diagrams for Cape Riva dacite determined by Cadoux et al. (2014) at 200 and 400 MPa to evaluate the role of vapour-saturated decompression in driving the growth of plagioclase rims and opx in the Therasia and Cape Riva magmas (Fig. 16). The phase diagram for the Cape Riva dacite is also applicable to the Therasia dacites, as their major element compositions are very similar. Cadoux et al. (2014) found that the Cape Riva dacite was stored prior to eruption at 875–900°C and vapour-saturated conditions with a molar fraction of H₂O (XH₂O) of 0.8–0.9 (Fig. 16a). Assuming that the dacite had ascended from a source at 400 MPa to pre-eruptive storage at 200 MPa at approximately constant XH₂O, this would trigger the growth of ~10% plagioclase and opx (also possibly a small amount of cpx), with plagioclase becoming stable before opx. These percentages are similar to the observed abundances of plagioclase and opx in the magmas (of which the plagioclase cores are volumetrically minor). The phase relationships are therefore consistent with decompression-driven growth of plagioclase and opx during ascent of dacitic melts from the mid-crustal storage region to the final level of pre-eruptive storage. Melt inclusions in plagioclase in Cape Riva pumices indeed record entrapment pressures ranging from 260 to 100 MPa (equivalent to depths of 10 to 4 km), which are consistent with phenocryst growth during ascent (Druitt et al. 2016). Cooling may also have played a part; however, it is difficult to cool large bodies

<figure><img src="image_20.png" /><figcaption>**Fig. 16.** Isobaric phase relationships of the Cape Riva dacite at (a) 200 MPa and (b) 400 MPa, at the FMQ oxygen buffer and as a function of the molar water fraction [XH₂O] of the coexisting fluid phase and of temperature (°C), redrawn from Cadoux et al. (2014). Curves are labelled with mineral names lying inside their stability fields. Plag, plagioclase; Opx, orthopyroxene; Cpx, clinopyroxene; A, amphibole; Ilm, ilmenite. Each black circle represents an experimental charge. Italic black numbers are charge crystallinities in weight per cent. Light grey dashed lines indicate the anorthite content (mol %) of the experimental plagioclase. Grey band shows pre-eruptive temperature range inferred from natural Fe–Ti oxides. (See text for explanation.)</figcaption></figure>

of magma over a few decades without transferring that
magma to a cooler environment.

Decompression without ascent, for example by out-
gassing of an over-pressurized magma reservoir, cannot
be ruled out as an alternative mechanism for driving
crystallization. However, the phase relationships show
that a pressure drop of 150–200 MPa is required to cause
~10 vol. % crystallization, and this would appear to be
too great for overpressure release at constant depth.
Overpressures in magma chambers typically do not ex-
ceed a few tens of MPa (Gudmundsson et al., 2012). The
decompression must therefore have involved significant
magma ascent. Plagioclase rims, unzoned plagioclase
crystals of rim composition, and opx crystals are envis-
aged to have nucleated and grown during this ascent.

We attribute the similarities in phenocryst textures,
compositions and residence timescales in the Therasia
and Cape Riva magmas to the ascent and decompression-
driven crystallization of each magma batch from the mid-
to upper crust, the difference being that the Therasia case
involved ascent of multiple, small volumes over 15 kyr,
whereas the Cape Riva case involved the ascent of much
larger volumes that amalgamated to form the final, pre-eruptive chamber. In each case rapid ascent-driven de-compression was followed by prolonged (years to decades) upper crustal storage prior to eruption.

This interpretation is consistent with the crystal resi-dence timescales. The plagioclase rims and opx in the Therasia dacites give NLRTs that are short compared with the mean repose time between Therasia eruptions. The mean time between eruptions is calculated as that between the youngest and oldest Therasia dacites, ~15 kyr, divided by 11, the minimum number of repose periods during that time (i.e. ~1.4 kyr per eruption). Even if we use the longer opx NLRT given by the diffu-sion coefficient of Dohmen et al. (2016), the longest NLRT are still shorter than the mean repose time (Fig. 15). The plagioclase rims and opx in the Therasia dacites therefore record growth and near-liquidus resi-dence over only a fraction of the mean inter-eruption re-pose time, and this is compatible with the rapid ascent, partial crystallization and eruption of multiple small vol-umes of dacite to build up the Therasia dome complex. In the case of the Cape Riva dacite, Fabbro et al. (2013) showed that assembly of the Cape Riva reservoir took place in less than 2800 ± 1400 years. Growth of the plagioclase rims and opx in the Cape Riva dacite is in-terpreted as possibly recording the last phases of as-sembly of this reservoir by the injection of new silicic melt, decades prior to eruption. The existence of two opx populations with different NLRT in the CR dacite (type c crystals have shorter NLRT than type a and b; Fig. 15) suggests that there may have been at least two pulses of reservoir assembly. Similarly short (decadal) phenocryst residence timescales have been estimated for other explosive eruptions (e.g. Morgan et al., 2006; Druitt et al., 2012; Saunders et al., 2012; Kilgour et al., 2014).

## Reconstruction of the events leading to the Cape Riva eruption

We combine our new petrological results with the field, chronological and geochemical results of Fabbro et al. (2013) to reconstruct the events during the transition from the interplinian (largely effusive) Therasia period to the Plinian (caldera-forming) Cape Riva eruption.

The period of silicic volcanism studied here followed the construction of the 350 m high (above present-day sea level) Skaros shield volcano, which occupied the present-day northern caldera basin. The shield was built by the repeated effusion of basaltic to andesitic lavas, and was complete by 54 ka. Following a period of low magma output, volcanism resumed at 39 ka with the dis-charge of ~2 km³ of mainly dacitic lavas to form the Therasia dome complex atop, and on the western flank of, Skaros. Fabbro et al. (2013) hypothesized that the lithostatic stress exerted by the Skaros edifice caused mantle-derived TH basalt to accumulate at depth (Pinel & Jaupart, 2000), possibly in the mid-crustal storage zone (~16 km) identified by Andújar et al. (2015), where it

generated TH-type silicic melts by a combination of frac-tional crystallization, crustal melting and defrosting of pre-existing plutons. These melts segregated from their parent crystal mushes to form crystal-poor silicic melts that ascended (along with entrained antecrystic plagio-clase) as a series of discrete volumes (Fig. 17). Some combination of cooling and vapour-saturated decompression led to the growth, during ascent, of opx crystals and of plagioclase rims on the pre-existing antecrystic cores. Successive melt volumes then either erupted to form the Therasia lavas or crystallized at depth as intrusions. Eruptions occurred within a few decades of ascent to the upper crustal holding level, as shown by the NLRT data. Small amounts of TH basaltic magma were entrained with some of the silicic melts and chilled to form enclaves; some also mixed with the dacite to form a hybrid andesite (Fabbro et al., 2013). A similar story of rapid ascent and crystallization of small magma volumes has been proposed at other volcanoes, such as Montserrat and the present-day Kameni edifice on Santorini (Zellmer et al., 1999, 2003). The compositional and thermal monotony of successive waves of Therasia melt over 15 kyr (Fabbro et al., 2013) may be due to their derivation from a mid-crustal storage region that was thermally buffered by a balance between the input of hot mafic magma, the release of latent heat through crystallization, and the loss of heat to the surrounding crust. The physical constraints of melt extraction may also have played a role (Dufek & Bachmann, 2010).

Sometime during the Therasia period a new batch of silicic melt, poorer in incompatible elements (CR dacite) began to ascend into the upper crust from a different region of mid-crustal intrusion (Fig. 17). The Fe/Mg ratio of the CR dacite is lower than those of the TH dacites, suggesting that it was sourced at slightly shallower levels in the mid-crust (Andújar et al., 2016). Most antecrystic plagioclase was derived from its own source mush, not that of the preceding Therasia melts, as demonstrated by their different contents of K and La (Fig. 9). The CR and TH dacites were probably stored in the same region of the shallow crust immediately before eruption, as the TH vents were located within the area that collapsed to form the Cape Riva caldera, and the similarity in the major-element compositions of the minerals found in the two dacites suggests that they were stored at comparable depths. The near-absence in Cape Riva pumice of cannibalized Therasia plagioclase suggests that any remaining Therasia-type dacite in the upper crust beneath the edifice must have largely crystallized prior to assembly of the upper crustal Cape Riva reservoir. The possible presence of the CR dacite as a chemical component in the 26 ka Cape Tripiti pumice layer of the Therasia complex (demonstrated by its relative depletion in incompatible trace elements compared with the rest of the TH dacites) shows that some CR dacite may have been residing in the upper crust by that time (Fabbro et al., 2013). However, the continued eruption of TH dacite from the same geographical footprint as the future Cape Riva caldera shows that most
<figure><img src="image_21.png" /><figcaption>**Fig. 17.** Schematic diagram summarizing the changes in the plumbing system of Santorini between 39 and 22 ka, showing the transition from dominantly effusive to caldera-forming states. Therasia-type basalt ponds in the middle crust, perhaps at a density interface, where it generates Therasia dacitic melts that ascend through the upper crust and erupt between 39 and 25 ka. Around 26 ka, the first Cape Riva type basalt is injected into the middle crust, where it begins to generate Cape Riva dacitic melt. Following the last Therasia eruption (~25 ka), runaway drainage of the Cape Riva mid-crustal body leads to assembly of the upper crustal Cape Riva magma chamber, which is then discharged explosively at $21.8 \pm 0.4$ ka. The model is based on the work on Fabbro *et al.* (2013), as well as the present paper. TH, Therasia; CR, Cape Riva.</figcaption></figure>

of the shallow Cape Riva reservoir must have been assembled over the $<2800 \pm 1400$ years following the last Therasia eruption. This implies drainage of $>10 \text{ km}^3$ of CR dacitic melt from its mid-crustal source volume at a time-averaged rate of $>0.004 \pm 0.002 \text{ km}^3 \text{ a}^{-1}$ between ~25 and $21.8 \pm 0.4$ ka (Fabbro *et al.*, 2013). Vapour-saturated growth of opx and plagioclase is consistent with phase equilibria (Fig. 16; Cadoux *et al.*, 2014) and is interpreted as having taken place during ascent into the upper crust at 4–6 km depth (Druitt *et al.*, 2016). Crystallization of opx with Al sector zoning is suggestive of rapid crystallization followed by residence in the upper crust prior to eruption. Injection of the new CR-type melts beneath the Skaros–Therasia edifice was sufficiently fast to allow runaway development of a melt-dominated magma reservoir (Gelman *et al.*, 2013; Schöpa & Annen, 2013). Although reservoir assembly may have taken many hundreds of years (but not much more than $2800 \pm 1400$ years), the short NLRT of the plagioclase margins and opx suggest that CR dacitic melt

was still being injected into the Cape Riva magma reservoir only years to decades prior to eruption.

## The transition from dominantly effusive to caldera-forming states at a caldera volcano

Better understanding of the processes governing transitions from interplinian to caldera-forming states at arc volcanoes is fundamental to hazard assessment. The traditional model of slow, incremental assembly of large magma reservoirs is being progressively challenged (Cashman & Sparks, 2013, and references therein). Seismic tomography studies at calderas commonly fail to detect large upper crustal melt volumes (Lees, 2007), suggesting that large upper crustal reservoirs are transient features on geological timescales. Petrological studies are providing a growing body of evidence for the assembly of caldera-forming magma reservoirs on timescales much shorter than the preceding interplinian period (Sutton *et al.*, 2000; Wilson & Charlier, 2009; Druitt *et al.*, 2012; Gualda *et al.*, 2012;
Solano et al., 2012; Allan et al., 2013, 2017; Chamberlain
et al., 2014). Owing to its good preservation, the
Therasia lava pile provides a window into the evolving
magmatic system leading up to the Cape Riva caldera-
forming eruption. The field, chronological and geo-
chemical data presented by Fabbro et al. (2013), along
with the petrological data of the present paper, favour
an interpretation in which the Cape Riva magma reser-
voir resulted from the rapid drainage of silicic melt from
a mid-crustal source region on a geologically short
timescale. The transition from interplinian to caldera-
forming conditions between 24-6 and 21-8 ka at
Santorini involved major reorganization of the crustal
plumbing system, with late-stage growth pulses of the
Cape Riva reservoir only decades to years prior to the
final eruption. This is similar to the events leading up to
the 3-6 ka Minoan eruption, as proposed by Druitt et al.
(2012). If the rapid growth of shallow silicic reservoirs is
common prior to caldera-forming eruptions, then geo-
physical and geochemical monitoring may warn us of
impending major eruptions at restless calderas.

ACKNOWLEDGEMENTS

We thank Jean-Luc Devidal for his expertise in electron
microprobe and ICP-MS analysis, and Jean-Marc Hénot
for his expertise in electron microscope imagery. Many
thanks go to the editor, Gerhard Wörner, and to re-
viewers Kathy Cashman, Tomiya Akihiko, Ralf
Gertisser, Michihiko Nakamura, Chiara Petrone and one
anonymous reviewer for their helpful comments.
Caroline Bouvet de Maisonneuve provided helpful com-
ments on an earlier version of this paper.

FUNDING

G.F. was supported by the French Région d'Auvergne.
F.C. was supported by the 'Magma Plumbing System'
project from Singapore's National Research
Foundation. The 3 weeks spent by G.F. at the Earth
Observatory of Singapore was supported by the
Merlion program of the Embassy of France in
Singapore. This is Laboratory of Excellence ClerVolc
contribution number 229.

SUPPLEMENTARY DATA

Supplementary data for this paper are available at
*Journal of Petrology* online.

REFERENCES

Allan, A. S. R., Morgan, D. J., Wilson, C. J. N. & Millet, M.-A. (2013). From mush to eruption in centuries: assembly of the super-sized Oruanui magma body. *Contributions to Mineralogy and Petrology* **166**, 143–164.

Allan, A. S. R., Barker, S. J., Millet, M.-A., Morgan, D. J., Rooyakkers, S. M., Schipper, C. I. & Wilson, C. J. N. (2017). A cascade of magmatic events during the assembly and

eruption of a super-sized magma body. *Contributions to Mineralogy and Petrology* **172**, 49.

Andújar, J., Scaillet, B., Pichavant, M. & Druitt, T. H. (2015). Differentiation conditions of a basaltic magma from Santorini, and its bearing on the production of andesite in arc settings. *Journal of Petrology* **56**, 765–794.

Andújar, J., Scaillet, B., Pichavant, M. & Druitt, T. H. (2016). Generation conditions of dacite and rhyodacite via the crystallization of an andesitic magma. Implications for the plumbing system at Santorini (Greece) and the origin of tholeiitic or calc-alkaline differentiation trends in arc magmas. *Journal of Petrology* **57**, 1887–1920.

Annen, C., Blundy, J. D. & Sparks, R. S. J. (2006). The genesis of intermediate and silicic magmas in deep crustal hot zones. *Journal of Petrology* **47**, 505–539.

Arce, J. L., Gardner, J. E. & Macías, J. L. (2013). Pre-eruptive conditions of dacitic magma erupted during the 21.7 ka Plinian event at Nevado de Toluca volcano, Central Mexico. *Journal of Volcanology and Geothermal Research* **249**, 49–65.

Bachmann, O. & Bergantz, G. (2008). The magma reservoirs that feed supereruptions. *Elements* **4**, 17–21.

Bacon, C. R. (1985). Implications of silicic vent patterns for the presence of large crustal magma chambers. *Journal of Geophysical Research* **90**, 11243–11252.

Bacon, C. R. & Druitt, T. H. (1988). Compositional evolution of the zoned calcalkaline magma chamber of Mount Mazama, Crater Lake, Oregon. *Contributions to Mineralogy and Petrology* **98**, 224–256.

Bacon, C. R. & Lowenstern, J. B. (2005). Late Pleistocene granodiorite source for recycled zircon and phenocrysts in rhyodacite lava at Crater Lake, Oregon. *Earth and Planetary Science Letters* **233**, 277–293.

Bailey, J. C., Jensen, E. S., Hansen, A., Kann, A. D. J. & Kann, K. (2009). Formation of heterogeneous magmatic series beneath North Santorini, South Aegean island arc. *Lithos* **110**, 20–36.

Bédard, J. H. (2006). Trace element partitioning in plagioclase feldspar. *Geochimica et Cosmochimica Acta* **70**, 3717–3742.

Bédard, J. H. (2007). Trace element partitioning coefficients between silicate melts and orthopyroxene: parameterizations of *D* variations. *Chemical Geology* **244**, 263–303.

Behrens, H., Johannes, W. & Schmalzried, H. (1990). On the mechanisms of cation diffusion processes in ternary feldspars. *Physics and Chemistry of Minerals* **17**, 62–78.

Bindeman, I. N., Davis, A. M. & Drake, M. J. (1998). Ion micro-probe study of plagioclase–basalt partition experiments at natural concentration levels of trace elements. *Geochimica et Cosmochimica Acta* **62**, 1175–1193.

Blundy, J. D. & Wood, B. J. (1991). Crystal-chemical controls on the partitioning of Sr and Ba between plagioclase feldspar, silicate melts, and hydrothermal solutions. *Geochimica et Cosmochimica Acta* **55**, 193–209.

Blundy, J. & Wood, B. (2003). Partitioning of trace elements between crystals and melts. *Earth and Planetary Science Letters* **210**, 383–397.

Boyce, J. W. & Hervig, R. L. (2008). Magmatic degassing histories from apatite volatile stratigraphy. *Geology* **36**, 63–66.

Brown, S. J. A. & Fletcher, I. R. (1999). SHRIMP U–Pb dating of the preeruption growth history of zircons from the 340 ka Whakamaru Ignimbrite, New Zealand: evidence for >250 k.y. magma residence times. *Geology* **27**, 1035–1038.

Cadoux, A., Scaillet, B., Druitt, T. H. & Deloule, E. (2014). Magma storage conditions of large Plinian eruptions of Santorini Volcano (Greece). *Journal of Petrology* **55**, 1129–1171.

Cashman, K. & Blundy, J. (2013). Petrological cannibalism: the chemical and textural consequences of incremental magma
body growth. *Contributions to Mineralogy and Petrology* **166**, 703–729.

Cashman, K. V. & Sparks, R. S. J. (2013). How volcanoes work: a 25 year perspective. *Geological Society of America Bulletin* **125**, 664–690.

Chamberlain, K. J., Wilson, C. J. N., Wooden, J. L., Charlier, B. L. A. & Ireland, T. R. (2014). New perspectives on the Bishop Tuff from zircon textures, ages and trace elements. *Journal of Petrology* **55**, 395–426.

Cherniak, D. J. (1995). Diffusion of lead in plagioclase and K-feldspar: an investigation using Rutherford Backscattering and Resonant Nuclear Reaction Analysis. *Contributions to Mineralogy and Petrology* **120**, 358–371.

Cherniak, D. J. (2002a). REE diffusion in feldspar. *Chemical Geology* **193**, 25–41.

Cherniak, D. J. (2002b). Ba diffusion in feldspar. *Geochimica et Cosmochimica Acta* **66**, 1641–1650.

Cherniak, D. J. (2003). Silicon self-diffusion in single-crystal natural quartz and feldspar. *Earth and Planetary Science Letters* **214**, 655–668.

Cherniak, D. J. (2010). Cation diffusion in feldspars. In: Zhang, Y. & Cherniak, D. J. (eds) *Diffusion in Minerals and Melts. Mineralogical Society of America and Geochemical Society, Reviews in Mineralogy and Geochemistry* **72**, 691–733.

Cherniak, D. J. & Dimanov, A. (2010). Diffusion in pyroxene, mica and amphibole. In: Zhang, Y. & Cherniak, D. J. (eds) *Diffusion in Minerals and Melts. Mineralogical Society of America and Geochemical Society, Reviews in Mineralogy and Geochemistry* **72**, 641–690.

Cherniak, D. J. & Watson, E. B. (1994). A study of strontium diffusion in plagioclase using Rutherford backscattering spectroscopy. *Geochimica et Cosmochimica Acta* **58**, 5179–5190.

Christensen, J. N. & DePaolo, D. J. (1993). Time scales of large volume silicic magma systems: Sr isotopic systematics of phenocrysts and glass from the Bishop Tuff, Long Valley, California. *Contributions to Mineralogy and Petrology* **113**, 100–114.

Cooper, K. M. & Kent, A. J. R. (2014). Rapid remobilization of magmatic crystals kept in cold storage. *Nature* **506**, 480–483.

Costa, F., Chakraborty, S. & Dohmen, R. (2003). Diffusion coupling between trace and major elements and a model for calculation of magma residence times using plagioclase. *Geochimica et Cosmochimica Acta* **67**, 2189–2200.

Costa, F., Dohmen, R. & Chakraborty, S. (2008). Time scales of magmatic processes from modelling the zoning patterns of crystals. In: Putirka, K. D. & Tepley, F. J., III (eds) *Minerals, Inclusions and Volcanic Processes. Mineralogical Society of America and Geochemical Society, Reviews in Mineralogy and Geochemistry* **69**, 545–594.

Crank, J. (1957). Two methods for the numerical solution of moving-boundary problems in diffusion and heat flow. *Quarterly Journal of Mechanics and Applied Mathematics* **10**, 220–231.

Davies, G. R., Halliday, A. N., Mahood, G. A. & Hall, C. M. (1994). Isotopic constraints on the production rates, crystallisation histories and residence times of pre-caldera silicic magmas, Long Valley, California. *Earth and Planetary Science Letters* **125**, 17–37.

Deer, W. A., Howie, R. A. & Zussman, J. (1992). *An Introduction to the Rock-Forming Minerals*. Pearson Education.

Dohmen, R., Ter Heege, J. H., Becker, H.-W. & Chakraborty, S. (2016). Fe–Mg interdiffusion in orthopyroxene. *American Mineralogist* **101**, 2210–2221.

Druitt, T. H. (1985). Vent evolution and lag breccia formation during the Cape Riva eruption of Santorini, Greece. *Journal of Geology* **93**, 439–454.

Druitt, T. H. & Bacon, C. R. (1989). Petrology of the zoned calc-alkaline magma chamber of Mount Mazama, Crater Lake, Oregon. *Contributions to Mineralogy and Petrology* **101**, 245–259.

Druitt, T. H., Edwards, L., Mellors, R. M., Pyle, D. M., Sparks, R. S. J., Lanphere, M., Davies, M. & Barriero, B. (1999). *Santorini Volcano. Geological Society, London, Memoirs* **19**.

Druitt, T. H., 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**, 77–80.

Druitt, T. H., Mercier, M., Florentin, L., Deloule, E., Cluzel, N., Médard, E. & Cadoux, A. (2016). Magma storage and extraction associated with Plinian and interplinian activity at Santorini Caldera (Greece). *Journal of Petrology* **57**, 461–494.

Dufek, J. & Bachmann, O. (2010). Quantum magmatism: magmatic compositional gaps generated by melt-crystal dynamics. *Geology* **38**, 687–690.

Endrun, B., Meier, T., Lebedev, S., Bohnhoff, M., Stavrakakis, G. & Harjes, H.-P. (2008). S velocity structure and radial anisotropy in the Aegean region from surface wave dispersion. *Geophysical Journal International* **174**, 593–616.

Faak, K., Chakraborty, S. & Coogan, L. A. (2013). Mg in plagioclase: experimental calibration of a new geothermometer and diffusion coefficients. *Geochimica et Cosmochimica Acta* **123**, 195–217.

Fabbro, G. N., Druitt, T. H. & Scaillet, S. (2013). Evolution of the crustal magma plumbing system during the build-up to the 22-ka caldera-forming eruption of Santorini (Greece). *Bulletin of Volcanology* **75**, 1–22.

Folkes, C. B., de Silva, S. L., Schmitt, A. K. & Cas, R. A. F. (2011). A reconnaissance of U–Pb zircon ages in the Cerro Galán system, NW Argentina: prolonged magma residence, crystal recycling, and crustal assimilation. *Journal of Volcanology and Geothermal Research* **206**, 136–147.

Frost, B. R. (1991). Introduction to oxygen fugacity and its petrologic importance. In: Lindsley, D. H. (ed.) *Oxide Minerals. Mineralogical Society of America, Reviews in Mineralogy* **25**, 1–9.

Ganguly, J. & Tazzoli, V. (1994). Fe²⁺–Mg interdiffusion in orthopyroxene; retrieval from the data on intracrystalline exchange reaction. *American Mineralogist* **79**, 930–937.

Gelman, S. E., Gutiérrez, F. J. & Bachmann, O. (2013). On the longevity of large upper crustal silicic magma reservoirs. *Geology* **41**, 759–762.

Gertisser, R., Preece, K. & Keller, J. (2009). The Plinian Lower Pumice 2 eruption, Santorini, Greece: magma evolution and volatile behaviour. *Journal of Volcanology and Geothermal Research* **186**, 387–406.

Giletti, B. J. & Casserly, J. E. D. (1994). Strontium diffusion kinetics in plagioclase feldspars. *Geochimica et Cosmochimica Acta* **58**, 3785–3793.

Giletti, B. J. & Shanahan, T. M. (1997). Alkali diffusion in plagioclase feldspar. *Chemical Geology* **139**, 3–20.

Grove, T. L., Baker, M. B. & Kinzler, R. J. (1984). Coupled CaAl-NaSi diffusion in plagioclase feldspar: experiments and applications to cooling rate speedometry. *Geochimica et Cosmochimica Acta* **48**, 2113–2121.

Gualda, G. A. R., Pamukcu, A. S., Ghiorso, M. S., Anderson, A. T., Sutton, S. R. & Rivers, M. L. (2012). Timescales of quartz crystallization and the longevity of the Bishop giant magma body. *PLoS One* **7**, e37492.

Gudmundsson, A., Kusumoto, S., Simmenes, T. H., Philipp, S. L., Larsen, B. & Lotveit, I. F. (2012). Effects of overpressure variations on fracture apertures and fluid transport. *Tectonophysics* **581**, 220–230.
Halliday, A. N., Mahood, G. A., Holden, P., Metz, J. M., Dempster, T. J. & Davidson, J. P. (1989). Evidence for long residence times of rhyolitic magma in the Long Valley magmatic system: the isotopic record in precaldera lavas of Glass Mountain. *Earth and Planetary Science Letters* **94**, 274–290.

Hildreth, W. (1979). The Bishop Tuff: evidence for the origin of compositional zonation in silicic magma chambers. In: Chapin, C. E. & Elston, W. E. (eds) *Ash-flow Tuffs*. Geological Society of America, Special Papers **180**, 43–76.

Hildreth, W. (1981). Gradients in silicic magma chambers: implications for lithospheric magmatism. *Journal of Geophysical Research: Solid Earth* **86**, 10153–10192.

Hildreth, W. & Wilson, C. J. N. (2007). Compositional zoning of the Bishop Tuff. *Journal of Petrology* **48**, 951–999.

Hollister, L. S. & Bence, A. E. (1967). Staurolite: sectoral compositional variations. *Science* **158**, 1053–1056.

Huber, C., Bachmann, O. & Manga, M. (2009). Homogenization processes in silicic magma chambers by stirring and mushification (latent heat buffering). *Earth and Planetary Science Letters* **283**, 38–47.

Huijsmans, J. P. P. (1985). Calc-alkaline lavas from the volcanic complex of Santorini, Aegean Sea, Greece. A petrological and stratigraphical study. PhD thesis, Rijksuniversiteit te Utrecht, Utrecht.

Huijsmans, J. P. P. & Barton, M. (1989). Polybaric geochemical evolution of two shield volcanoes from Santorini, Aegean Sea, Greece: evidence for zoned magma chambers from cyclic compositional variations. *Journal of Petrology* **30**, 583–625.

Huijsmans, J. P. P., Barton, M. & Salters, V. J. M. (1988). Geochemistry and evolution of the calc-alkaline volcanic complex of Santorini, Aegean Sea, Greece. *Journal of Volcanology and Geothermal Research* **34**, 283–306.

Jolivet, L., Faccenna, C., Huet, B., Labrousse, L., Le Pourhiet, L., Lacombe, O., Lecomte, E., Burov, E., Denèle, Y., Brun, J.-P., Philippon, M., Paul, A., Salaün, G., Karabulut, H., Piromallo, C., Monié, P., Gueydan, F., Okay, A. I., Oberhäsli, R., Pourteau, A., Augier, R., Gadenne, L. & Driussi, O. (2013). Aegean tectonics: strain localisation, slab tearing and trench retreat. *Tectonophysics* **597–598**, 1–33.

Karagianni, E. E., Papazachos, C. B., Panagiotopoulos, D. G., Suhadolc, P., Vuan, A. & Panza, G. F. (2005). Shear velocity structure in the Aegean area obtained by inversion of Rayleigh waves. *Geophysical Journal International* **160**, 127–143.

Kilgour, G. N., Saunders, K. E., Blundy, J. D., Cashman, K. V., Scott, B. J. & Miller, C. A. (2014). Timescales of magmatic processes at Ruapehu volcano from diffusion chronometry and their comparison to monitoring data. *Journal of Volcanology and Geothermal Research* **288**, 62–75.

Klügel, A. (2001). Prolonged reactions between harzburgite xenoliths and silica-undersaturated melt: implications for dissolution and Fe–Mg interdiffusion rates of orthopyroxene. *Contributions to Mineralogy and Petrology* **141**, 1–14.

Konstantinou, K. I. (2010). Crustal rheology of the Santorini–Amorgos zone: implications for the nucleation depth and rupture extent of the 9 July 1956 Amorgos earthquake, southern Aegean. *Journal of Geodynamics* **50**, 400–409.

Koyaguchi, T. & Kaneko, K. (1999). A two-stage thermal evolution model of magmas in continental crust. *Journal of Petrology* **40**, 241–254.

Koyaguchi, T. & Kaneko, K. (2000). Thermal evolution of silicic magma chambers after basalt replenishments. *Transactions of the Royal Society of Edinburgh: Earth Sciences* **91**, 47–60.

Krimer, D. & Costa, F. (2017). Evaluation of the effects of 3D diffusion, crystal geometry, and initial conditions on retrieved time-scales from Fe–Mg zoning in natural oriented orthopyroxene crystals. *Geochimica et Cosmochimica Acta* **196**, 271–288.

Kudo, A. M. & Weill, D. F. (1970). An igneous plagioclase thermometer. *Contributions to Mineralogy and Petrology* **25**, 52–65.

Lasaga, A. C. (1998). *Kinetic Theory in the Earth Sciences*. Princeton University Press.

LaTourrette, T. & Wasserburg, G. J. (1998). Mg diffusion in anorthite: implications for the formation of early solar system planetesimals. *Earth and Planetary Science Letters* **158**, 91–108.

Lees, J. M. (2007). Seismic tomography of magmatic systems. *Journal of Volcanology and Geothermal Research* **167**, 37–56.

Le Pichon, X. & Angelier, J. (1979). The Hellenic arc and trench system: a key to the neotectonic evolution of the eastern Mediterranean area. *Tectonophysics* **60**, 1–42.

Liu, M. & Yund, R. A. (1992). NaSi-CaAl interdiffusion in plagioclase. *American Mineralogist* **77**, 275–283.

Liu, Y., Anderson, A. T., Wilson, C. J. N., Davis, A. M. & Steele, I. M. (2005). Mixing and differentiation in the Oruanui rhyolitic magma, Taupo, New Zealand: evidence from volatiles and trace elements in melt inclusions. *Contributions to Mineralogy and Petrology* **151**, 71–87.

Longpré, M.-A., Stix, J., Costa, F., Espinoza, E. & Muñoz, A. (2014). Magmatic processes and associated timescales leading to the January 1835 eruption of Cosigüina Volcano, Nicaragua. *Journal of Petrology* **55**, 1173–1201.

Mahood, G. A. (1990). Second reply to comment of R.S.J. Sparks, H.E. Huppert and C.J.N. Wilson on 'Evidence for long residence times of rhyolitic magma in the Long Valley magmatic system: the isotopic record in the precaldera lavas of Glass Mountain'. *Earth and Planetary Science Letters* **99**, 395–399.

Martel, C. & Schmidt, B. C. (2003). Decompression experiments as an insight into ascent rates of silicic magmas. *Contributions to Mineralogy and Petrology* **144**, 397–415.

Metz, J. M. & Mahood, G. A. (1985). Precursors to the Bishop Tuff Eruption: Glass Mountain, Long Valley, California. *Journal of Geophysical Research* **90**, 11121–11126.

Metz, J. M. & Mahood, G. A. (1991). Development of the Long Valley, California, magma chamber recorded in precaldera rhyolite lavas of Glass Mountain. *Contributions to Mineralogy and Petrology* **106**, 379–397.

Michaud, V., Clocchiatti, R. & Sbrana, S. (2000). The Minoan and post-Minoan eruptions, Santorini (Greece), in the light of melt inclusions: chlorine and sulphur behaviour. *Journal of Volcanology and Geothermal Research* **99**, 195–214.

Morgan, D. J., Blake, S., Rogers, N. W., De Vivo, B., Rolandi, G. & Davidson, J. P. (2006). Magma chamber recharge at Vesuvius in the century prior to the eruption of A.D. 79. *Geology* **34**, 845–848.

Nakagawa, K., Nagahara, H., Ozawa, K., Tachibana, S. & Yasuda, A. (2005). Experimental determination of AlAI-SiMg interdiffusion coefficient in orthopyroxene. *Japan Earth and Planetary Science Joint Meeting, Tokyo*, 22–26 May 2005, abstract K038-013.

Petrone, C. M., Bugatti, G., Braschi, E. & Tommasini, S. (2016). Pre-eruptive magmatic processes re-timed using a non-isothermal approach to magma chamber dynamics. *Nature Communications* **7**, 12946.

Pinel, V. & Jaupart, C. (2000). The effect of edifice load on magma ascent beneath a volcano. *Philosophical
*Transactions of the Royal Society of London, Series A* **358**, 1515–1532.

Putirka, K. D. (2005). Igneous thermometers and barometers based on plagioclase + liquid equilibria: tests of some existing models and new calibrations. *American Mineralogist* **90**, 336–346.

Rasband, W. S. (1997). *ImageJ*. US National Institutes of Health.

Ruprecht, P. & Cooper, K. M. (2012). Integrating the uranium-series and elemental diffusion geochronometers in mixed magmas from Volcán Quizapu, central Chile. *Journal of Petrology* **53**, 841–871.

Sachpazi, M., Hirn, A., Nercessian, A., Avedik, F., McBride, J., Loucoyannakis, M., Nicolich, R. & the STREAMERS-PROFILES group (1997). A first coincident normal-incidence and wide-angle approach to studying the extending Aegean crust. *Tectonophysics* **270**, 301–312.

Saunders, K., Blundy, J., Dohmen, R. & Cashman, K. (2012). Linking petrology and seismology at an active volcano. *Science* **336**, 1023–1027.

Saunders, K. E., Morgan, D. J., Baker, J. A. & Wysoczanski, R. J. (2010). The magmatic evolution of the Whakamaru supereruption, New Zealand, constrained by a microanalytical study of plagioclase and quartz. *Journal of Petrology* **51**, 2465–2488.

Scaillet, B. & Evans, B. W. (1999). The 15 June 1991 eruption of Mount Pinatubo. I. Phase equilibria and pre-eruption P–T–fO$_2$–fH$_2$O conditions of the dacite magma. *Journal of Petrology* **40**, 381–411.

Schöpa, A. & Annen, C. (2013). The effects of magma flux variations on the formation and lifetime of large silicic magma chambers. *Journal of Geophysical Research: Solid Earth* **118**, 926–942.

Schwandt, C. S. & McKay, G. A. (2006). Minor- and trace-element sector zoning in synthetic enstatite. *American Mineralogist* **91**, 1607–1615.

Schwandt, C. S., Cygan, R. T. & Westrich, H. R. (1998). Magnesium self-diffusion in orthoenstatite. *Contributions to Mineralogy and Petrology* **130**, 390–396.

Schwarz, M. (2000). Tephrakorrelation im östlichen Mittelmeer (Meteor M40/4 Kerne). Diploma Thesis, Albert-Ludwigs-Universität Freiburg i. Br, Freiburg.

Self, S. & Blake, S. (2008). Consequences of explosive supereruptions. *Elements* **4**, 41–46.

Sigurdsson, H., Carey, S. & Devine, J. D. (1990). Assessment of mass, dynamics and environmental effects of the Minoan eruption of Santorini Volcano. In: Hardy, D. A., Keller, J., Galanopoulos, V. P., Flemming, N. C. & Druitt, T. H. (eds) *Thera and the Aegean World III*. Thera Foundation, pp. 100–112.

Simon, J. I. & Reid, M. R. (2005). The pace of rhyolite differentiation and storage in an 'archetypical' silicic magma system, Long Valley, California. *Earth and Planetary Science Letters* **235**, 123–140.

Smith, D. & Barron, B. R. (1991). Pyroxene-garnet equilibration during cooling in the mantle. *American Mineralogist* **76**, 1950–1963.

Solano, J. M. S., Jackson, M. D., Sparks, R. S. J., Blundy, J. D. & Annen, C. (2012). Melt segregation in deep crustal hot zones: a mechanism for chemical differentiation, crustal assimilation and the formation of evolved magmas. *Journal of Petrology* **53**, 1999–2026.

Sparks, R. S. J., Huppert, H. E. & Wilson, C. J. N. (1990). Comment on 'Evidence for long residence times of rhyolitic magma in the Long Valley magmatic system: the isotopic record in precaldera lavas of Glass Mountain' by A.N. Halliday, G.A. Mahood, P. Holden, J.M. Metz, T.J. Dempster and J.P. Davidson. *Earth and Planetary Science Letters* **99**, 387–389.

Stix, J., Goff, F., Gorton, M. P., Heiken, G. & Garcia, S. R. (1988). Restoration of compositional zonation in the Bandelier silicic magma chamber between two caldera-forming eruptions: geochemistry and origin of the Cerro Toledo Rhyolite, Jemez Mountains, New Mexico. *Journal of Geophysical Research: Solid Earth* **93**, 6129–6147.

Sutton, A. N., Blake, S. & Wilson, C. J. N. (1995). An outline geochemistry of rhyolite eruptives from Taupo volcanic centre, New Zealand. *Journal of Volcanology and Geothermal Research* **68**, 153–175.

Sutton, A. N., Blake, S., Wilson, C. J. N. & Charlier, B. L. A. (2000). Late Quaternary evolution of a hyperactive rhyolite magmatic system: Taupo volcanic centre, New Zealand. *Journal of the Geological Society, London* **157**, 537–552.

Tirel, C., Gueydan, F., Tiberi, C. & Brun, J.-P. (2004). Aegean crustal thickness inferred from gravity inversion. Geodynamical implications. *Earth and Planetary Science Letters* **228**, 267–280.

Tomiya, A. & Takahashi, E. (2005). Evolution of the magma chamber beneath Usu Volcano since 1663: a natural laboratory for observing changing phenocryst compositions and textures. *Journal of Petrology* **46**, 2395–2426.

Vaggelli, G., Pellegrini, M., Vougioukalakis, G., Innocenti, S. & Francalanci, L. (2009). Highly Sr radiogenic tholeiitic magmas in the latest inter-Plinian activity of Santorini volcano, Greece. *Journal of Geophysical Research* **114**, B06201.

Van Orman, J. A., Cherniak, D. J. & Kita, N. T. (2014). Magnesium diffusion in plagioclase: Dependence on composition, and implications for thermal resetting of the $^{26}$Al-$^{26}$Mg early solar system chronometer. *Earth and Planetary Science Letters* **385**, 79–88.

Vazquez, J. A. & Reid, M. R. (2004). Probing the accumulation history of the voluminous Toba magma. *Science* **305**, 991–994.

Wark, D. A., Hildreth, W., Spear, F. S., Cherniak, D. J. & Watson, E. B. (2007). Pre-eruption recharge of the Bishop magma system. *Geology* **35**, 235–238.

Welsch, B., Faure, F., Famin, V., Baronnet, A. & Bachelery, P. (2013). Dendritic crystallization: a single process for all the textures of olivine in basalts? *Journal of Petrology* **54**, 539–574.

Wilcock, J., Goff, F., Minarik, W. G. & Stix, J. (2012). Magmatic recharge during the formation and resurgence of the Valles Caldera, New Mexico, USA: evidence from quartz compositional zoning and geothermometry. *Journal of Petrology* **53**, 635–664.

Wilson, C. J. N. & Charlier, B. L. A. (2009). Rapid rates of magma generation at contemporaneous magma systems, Taupo Volcano, New Zealand: insights from U–Th model-age spectra in zircons. *Journal of Petrology* **50**, 875–907.

Wotzlaw, J.-F., Schaltegger, U., Frick, D. A., Dungan, M. A., Gerdes, A. & Günther, D. (2013). Tracking the evolution of large-volume silicic magma reservoirs from assembly to supereruption. *Geology* **41**, 867–870.

Zellmer, G. F., Blake, S., Vance, D., Hawkesworth, C. & Turner, S. (1999). Plagioclase residence times at two island arc volcanoes (Kameni Islands, Santorini, and Soufriere, St. Vincent) determined by Sr diffusion systematics. *Contributions to Mineralogy and Petrology* **136**, 345–357.

Zellmer, G. F., Sparks, R. S. J., Hawkesworth, C. J. & Wiedenbeck, M. (2003). Magma emplacement and remobilization timescales beneath Montserrat: insights from Sr and Ba zonation in plagioclase phenocrysts. *Journal of Petrology* **44**, 1413–1431.
# APPENDIX: DIFFUSION MODELLING

## MODELLING MG DIFFUSION IN PLAGIOCLASE

### Initial conditions

To model the diffusion of Mg in plagioclase, it is first necessary to estimate the initial Mg concentration profiles of the crystals. This was done in two steps. First, an equation relating the Ti and Mg concentrations in Santorini whole-rock data was constructed; this was assumed to approximate the liquid line of descent, as confirmed by recent melt inclusion studies (Druitt et al., 2016). Then, a similar relationship between Ti and An was found, and the two equations were combined to create the final equation for the initial Mg content in terms of An. This assumes that Ti does not diffuse significantly, so that the Ti content of the plagioclase represents that of the original melt. This assumption is confirmed by our results.

The relationship between Ti and Mg in the whole-rock data can be fitted by two exponential functions: one for the rocks produced by Fe-Ti oxide-absent fractional crystallization (where Ti increases with the degree of fractionation), and one for where Fe-Ti oxides are present during fractionation (where Ti decreases with the degree of fractionation; Fig. 18). Both of these equations have the form

$$C_{\text{Mg}}^{\text{liq}} = a_1 \exp(a_2 C_{\text{Ti}}^{\text{liq}}) \quad (3)$$

where $a_1$ and $a_2$ are coefficients to be determined. These can be combined with linear functions for Ti against the molar anorthite fraction, $X_{An}$, where Ti

concentrations are those in a melt in equilibrium with the plagioclase:

$$C_{\text{Ti}}^{\text{liq}} = b_1 X_{\text{An}} + b_2 \quad (4)$$

where $b_1$ and $b_2$ are also coefficients to be determined.

Combining these two equations, we obtain the concentration of Mg in the melts that crystallized plagioclase of a particular anorthite content:

$$C_{\text{Mg}}^{\text{liq}} = a_1 \exp[a_2(b_1 X_{\text{An}} + b_2)]. \quad (5)$$

Using the partition coefficients of Bindeman et al. (1998), we can calculate the concentration of Mg in the plagioclase that would be in equilibrium with these melts:

$$\begin{align}
C_{\text{Mg}}^{\text{xl}} &= a_1 \exp(a_2 b_1 X_{\text{An}} + a_2 b_2) \exp\left(\frac{AX_{\text{An}} + B}{RT}\right) \tag{6} \\
&= a_1 \exp\left[a_2 b_2 + \left(a_2 b_1 + \frac{A}{RT}\right) X_{\text{An}} + \frac{B}{RT}\right]. \nonumber
\end{align}$$

Substituting the numbers into this equation does not quite fit the data: the calculated initial concentrations are lower than those measured in the rims of the plagioclase. This may be a problem with the partition coefficients: Bindeman et al. (1998) measured theirs in basalt whereas our plagioclase grew from a dacite, and melt composition has been shown to affect trace element partitioning (Blundy & Wood, 2003). The same issue has also been noted by Druitt et al. (2012) and Longpré et al. (2014). Multiplying the equation by a constant fixes this, however. The final equation used to predict the initial concentration of Mg in plagioclase is (Fig. 19)

<figure><img src="image_22.png" /><figcaption>Fig. 18. (a) Whole-rock (filled symbols), groundmass (open symbols) and enclave (stars) compositions of the Therasia, Cape Riva, Skaros and Upper Scoria 2 (US2) eruptions (Huijsmans, 1985; Druitt et al., 1999; Fabbro et al., 2013). (b) Melt equivalent compositions for plagioclases from the Cape Riva and Therasia eruptions. Continuous black line is the equation used to estimate the initial Mg concentration of the plagioclase; the grey band is an estimate of the uncertainty covering about two-thirds of the data. Each blue field in (b) contains all the analyses from the profiles in an individual crystal in the Therasia dacites, and the green fields are individual crystals from the Cape Riva dacite.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Fig. 19. Estimate of the initial Mg concentration of plagioclase crystals, shown by the continuous black line, from equation (7). The uncertainty, from equations 8 and 9, is shown by the grey band. Each blue field contains all the analyses from the profiles in an individual crystal in the Therasia dacites, and the green fields are individual crystals from the Cape Riva dacite. Symbols are as in Figure 18.</figcaption></figure>

$$C_{\text{Mg}}^{\text{init}} = \begin{cases} 310 \exp \left[ \left( 6.13 + \frac{A}{RT} \right) X_{\text{An}} + \frac{B}{RT} \right] & \text{for An} < 69 \\ 2130 \exp \left[ \left( 3.33 + \frac{A}{RT} \right) X_{\text{An}} + \frac{B}{RT} \right] & \text{for An} \ge 69. \end{cases} \quad (7)$$

The potential range of initial Mg concentrations is given by the following bounds, which encompass about two-thirds of the data:

$$C_{\text{Mg}}^{\text{init-low}} = \begin{cases} 340 \exp \left[ \left( 5.62 + \frac{A}{RT} \right) X_{\text{An}} + \frac{B}{RT} \right] & \text{for An} < 72 \\ 2130 \exp \left[ \left( 3.33 + \frac{A}{RT} \right) X_{\text{An}} + \frac{B}{RT} \right] & \text{for An} \ge 72 \end{cases} \quad (8)$$

$$C_{\text{Mg}}^{\text{init-high}} = \begin{cases} 276 \exp \left[ \left( 6.72 + \frac{A}{RT} \right) X_{\text{An}} + \frac{B}{RT} \right] & \text{for An} < 66 \\ 3150 \exp \left[ \left( 3.03 + \frac{A}{RT} \right) X_{\text{An}} + \frac{B}{RT} \right] & \text{for An} \ge 66. \end{cases} \quad (9)$$

## Equilibrium profiles

Equilibrium profiles for the plagioclase crystals were calculated using the plagioclase–melt partition function of Bindeman et al. (1998). As the plagioclase–melt partition coefficient of element *i* is defined as

$$K_i^{\text{plag/liq}} = \frac{C_i^{\text{plag}}}{C_i^{\text{liq}}} \quad (10)$$

we can combine the partition coefficients for two different plagioclase compositions (plag-1 and plag-2) as follows:

$$\frac{K_i^{\text{plag-1/liq}}}{K_i^{\text{plag-2/liq}}} = \frac{C_i^{\text{plag-1}}/C_i^{\text{liq}}}{C_i^{\text{plag-2}}/C_i^{\text{liq}}} = \frac{C_i^{\text{plag-1}}}{C_i^{\text{plag-2}}} = K_i^{\text{plag-1-plag-2}}. \quad (11)$$

If the rim is assumed to be in equilibrium with the melt, then the equilibrium profile can be constructed by working inwards. We can use equation (11) to calculate the concentration of *i* in each part of the crystal that is in equilibrium with the part of the crystal immediately rimward of it. Diffusion will tend to drive the composition of the plagioclase towards this equilibrium profile (Costa et al., 2003).

## Diffusion coefficient of Mg in plagioclase

Three studies have systematically investigated the diffusion coefficient of Mg in plagioclase. LaTourrette & Wasserburg (1998) carried out self-diffusion experiments in crystals of An₉₅ at 1200–1400°C, Faak et al. (2013) carried out chemical diffusion experiments with plagioclase crystals of An₅₀ to An₆₇ at temperatures of 1100–1200°C, and Van Orman et al. (2014) carried out chemical diffusion experiments with plagioclase crystals with compositions between An₂₃ and An₉₃ at temperatures of 800–1150°C. All three studies reported similar activation energies: LaTourrette & Wasserburg (1998) gave 254 ± 43 and 278 ± 43 kJ mol⁻¹ for diffusion parallel to the *b*-axis and *c*-axis, respectively; Faak et al. (2013) gave 321 kJ mol⁻¹; and Van Orman et al. (2014) gave 287 ± 10 kJ mol⁻¹. Both LaTourrette & Wasserburg (1998) and Van Orman et al. (2014) found slight anisotropy for Mg diffusion, with diffusion parallel to the *c*-axis up to 2–3 times faster than that parallel to the *b*-axis. Van Orman et al. (2014) found a large compositional (i.e. Xₐₙ) dependence, similar to that reported for Sr, Pb, Ba, Nd and Ca (Behrens et al., 1990; Cherniak & Watson, 1994; Cherniak, 1995, 2002a, 2002b; LaTourrette & Wasserburg, 1998). However, Faak et al. (2013) did not find this compositional dependence over the range of plagioclase compositions that they studied (An₅₀–₆₇). Faak et al. (2013) reported a dependence on the activity of silica, with Dₘ₅plag ∝ (a₅SiO₅)²·⁶.

We use the data of Van Orman et al. (2014), as their experiments cover the largest range of temperatures and plagioclase compositions. Importantly, their temperature range includes the temperatures of Santorini dacites and rhyodacites (850–930°C; Sigurdsson et al., 1990; Fabbro et al., 2013; Cadoux et al., 2014), so that extrapolation to lower temperatures is not required. Their complete equation for the diffusion coefficient of Mg in plagioclase, Dₘ₅plag (m² s⁻¹), as a function of temperature and plagioclase composition is

$$D_{\text{Mg}}^{\text{plag}} = \exp \left[ \left( -6.06 \pm 1.10 \right) - \left( 7.96 \pm 0.42 \right) X_{\text{An}} - \frac{287000 \pm 10000}{RT} \right]. \quad (12)$$

## Numerical modelling

The diffusion of Mg in plagioclase was modelled numerically using the method of Costa et al. (2003). The changes in concentration owing to diffusion of a trace
element in one dimension in plagioclase such as Mg, C<sub>Mg</sub>, are given by

$$ \frac{\partial C_{Mg}}{\partial t} = \frac{\partial}{\partial x} \left( D_{Mg} \frac{\partial C_{Mg}}{\partial x} - D_{Mg} C_{Mg} \frac{A}{RT} \frac{\partial X_{An}}{\partial x} \right) \quad (13) $$

where D<sub>Mg</sub> is the diffusion coefficient of Mg in plagioclase, and A is the parameter that describes the variation of the partition coefficient of Mg with An content in equation (1). Equation (13) was then expanded and converted to a finite-difference scheme, and was implemented using Matlab scripts. The rimward end of the profile was fixed assuming the first measured point was in equilibrium with the melt it resided in. The profile was assumed to be symmetrical by implicitly mirroring it around the innermost point, except where the profile was significantly longer than half the crystal. Where the profile spanned the entire crystal, both ends of the profile were assumed to be in equilibrium with the melt.

The diffusion times of all the profiles were first estimated using one-stage models: the whole crystal was assumed to have grown instantaneously, with no diffusion, then allowed to re-equilibrate with the melt at 880°C (the typical Fe-Ti oxide temperature of the magmas concerned; Fabbro et al., 2013; Cadoux et al., 2014). The time taken for the calculated profile to match the measured profile was then assumed to be the residence time of the crystal at magmatic temperatures. Some of the profiles could, however, be better fitted by a two-stage model. The two-stage models were run using the technique of Druitt et al. (2012). The crystal was divided into two visually distinct parts based on its zoning pattern. (1) The core of the crystal was assumed to grow instantaneously, then allowed to partially re-equilibrate. The outermost boundary was assumed to still have its initial Mg concentration, and to be in equilibrium with its host melt. Diffusion was modelled at temperatures between 880 and 930°C, depending on the crystal An content [equation (2)]. (2) The rim of the crystal was then assumed to grow instantaneously, and the whole crystal was allowed to re-equilibrate diffusively with the melt, as in the single-stage models. Although this approach is still a simplification of the processes of concurrent growth and diffusion, it is a closer approximation to reality than one-stage models. This method allows two timescales to be calculated: the residence time of the core of the crystal before the growth of the rim, and the residence time of the whole crystal, both at the temperatures appropriate to the respective An contents.

## Uncertainties

Uncertainties in the near-liquidus residence times owing to data fitting were quantitatively assessed as follows. As the models were run, the time at which the calculated profile first reached the measured Mg concentrations was noted. The diffusion model was

then allowed to continue until the calculated profile fell below most of the measured Mg concentrations. Although this approach is somewhat subjective, it was conservatively applied to produce the largest possible range of ages supported by the data.

Another major source of uncertainty in the diffusion models is the initial conditions. Three profiles were selected to investigate the effect of the initial Mg concentration on the residence times recovered from the diffusion models. The three profiles were run starting at the three different estimates of the initial concentration ['best', 'low', and 'high', as given by equations (7), (8), and (9), respectively]. The time each profile took to reach different levels of equilibration was then compared for the different starting conditions. The effect of using either the high or low initial Mg concentrations, rather than the 'best' initial concentrations, was to change the estimated near-liquidus residence times by between 15 and 65 years (with an average difference of 40 years).

We can also assess the impact that the uncertainty in the initial conditions has on the near-liquidus residence times by looking at the time it takes for the calculated profile to drop from the 'best' initial profile to below the low bound of the initial concentration. In most cases this is less than 40 years. The exceptions are in the centres of long plateaux of the Mg concentration and regions affected by 'uphill diffusion'. For example, the very centre of GS11-30b XL15 (Supplementary Data Fig. S1.10) stays at roughly the same Mg concentration for the first 180 years. This is due to the fact that the diffusion front has to migrate to the centre. In the case of 'uphill diffusion' (e.g. GS10-14 XL32, Supplementary Data Fig. S1.3), it is Mg diffusing out of regions with high initial Mg concentrations that stops the Mg profile from dropping (even driving the Mg concentrations 'uphill', away from equilibrium). In both cases, changing the initial Mg concentration does not change this behaviour; instead it adds or subtracts ~40 years on top of the time it takes for the Mg concentrations in the regions affected to start dropping. A value of ±40 years was taken as the uncertainty owing to potential variations in the initial Mg concentration, and was added to the range of times estimated from the range of possible fits to the data (as described above). In crystals with near-liquidus residence times shorter than 40 years, however, the time taken for the one-stage models to drop from the best estimate of the initial concentration to below the lower bound of the estimates of the initial Mg concentrations was used instead.

The final source of uncertainty in the models is the uncertainty in the diffusion coefficient. Because time is inversely proportional to the diffusion coefficient, we can assess this analytically. The biggest source of error in the diffusion coefficient is the temperature at which the diffusion is assumed to have taken place. The
relationship between the diffusion coefficient and the temperature takes the form

$$D = D_0 \exp\left(\frac{-Q}{RT}\right). \qquad (14)$$

The ratio between the time diffusion takes at two different temperatures, T₁ and T₂, can therefore be written as

$$\frac{t_1}{t_2} = \frac{D_0 \exp\left(-\frac{Q}{RT_2}\right)}{D_0 \exp\left(-\frac{Q}{RT_1}\right)} \qquad (15)$$

$$\ln\left(\frac{t_1}{t_2}\right) = \frac{Q}{RT_1} - \frac{Q}{RT_2}. \qquad (16)$$

Fe–Ti oxides from two of the three analysed pumice fall deposits from the Therasia dome complex, along with the Cape Riva eruption, give temperatures of within about ±25°C of 880°C (Fabbro et al., 2013; Cadoux et al., 2014). This temperature range is therefore used to assess the uncertainty in the diffusion models.

Using equation (12) for An₄₀, diffusion at 855°C is ~2 times as long as at 880°C, and diffusion at 905°C is ~2 times faster. This approach may not be strictly applicable, however, as the diffusion coefficient of Mg in plagioclase is not constant along the length of the crystal. This leads to a ∂D/∂x term in the expansion of the diffusion equation. However, running models at 855 and 905°C confirms that a factor of two is a reasonable estimate of the uncertainty in the models owing to variation in temperature. The lower and upper bounds on the range of near-liquidus residence times estimated from the possible fits to the data and the uncertainties in the initial conditions were therefore multiplied or divided by a factor of two to take this into account.

## MODELLING Fe–Mg DIFFUSION IN ORTHOPYROXENE

### Initial conditions

The initial zone boundaries were assumed to be sharp step functions, which gives the maximum times that the observed Mg# variations could survive at high temperature. The initial Mg# on either side of the boundary was estimated by taking the values of Mg# at a distance away from the boundary where they appeared to be unmodified by diffusion. The weakly zoned orthopyroxene (opx) crystals do not have obvious zone boundaries; however, many of them do have small variations in Mg#. These variations often have profile shapes that are typical of diffusion between two distinct zones.

### Diffusion coefficient of Fe–Mg exchange in orthopyroxene

Fe–Mg interdiffusion in opx varies with opx composition, oxygen fugacity and temperature (Ganguly & Tazzoli, 1994; Klügel, 2001). The equation for the diffusion coefficient, D<sub>Fe–Mg</sub>, used for the models presented here is the

same as that used by Allan et al. (2013), corrected for a typographic error and with D in m² s⁻¹:

$$D_{\text{Fe-Mg}}^{\text{Opx}} = \left[ \frac{fO_2(\text{sample})}{fO_2(\text{IW buffer})} \right]^{1/6} \times 10^{-[9.54+2.6(1-X_{\text{Mg}})-(12530/T)]} \qquad (17)$$

where X<sub>Mg</sub> is the molar Mg fraction [Mg/(Mg + Fe)], and T is the temperature in K. For comparison, we also show the results using the equation of Dohmen et al. (2016):

$$D_{\text{Fe-Mg}}^{\text{Opx}} = 1.12 \times 10^{-6} \times (fO_2)^{0.053} \times \exp\left(\frac{-308000}{RT}\right) \times 10^{X_{\text{Fe}}-0.09} \qquad (18)$$

where R is the molar gas constant. Fe–Ti oxides from the Cape Riva dacite and two of the three Therasia pumice fall deposits analysed plot along the fayalite–magnetite–quartz (FMQ) oxygen buffer (Fabbro et al., 2013; Cadoux et al., 2014), which is about four log units higher than the iron–wüstite (IW) buffer used in the calibration of the equation above (Frost, 1991). The diffusion coefficient is strongly dependent on the temperature, which for the Cape Riva dacite and two of the three analysed Therasia dacites is about 880 ± 25°C (Fabbro et al., 2013; Cadoux et al., 2014). The models were therefore run at 880°C, and at the FMQ oxygen buffer. The effect of using different temperatures and oxygen fugacities is discussed below.

### Numerical modelling

Because the variation in X<sub>Mg</sub> in each profile was always small (<0.1), the dependence of the Mg–Fe diffusion coefficient on X<sub>Mg</sub> within a single profile was ignored. This means that the profiles can be fitted by a simple analytical solution to the diffusion equation (Crank, 1957):

$$X_{\text{Mg}}(x, t) = X_{\text{Mg}}^{\text{left}} + \frac{X_{\text{Mg}}^{\text{left}} + X_{\text{Mg}}^{\text{right}}}{2} \left[ 1 + \text{erf}\left(\frac{x}{2\sqrt{Dt}}\right) \right] \qquad (19)$$

where X<sub>Mg</sub><sup>left</sup> and X<sub>Mg</sub><sup>right</sup> are the initial molar Mg fraction to the left and the right of the step function, respectively; x is the distance, centred on the initial step function; t is the time; and erf(u) is the error function. The diffusion coefficient was calculated using an X<sub>Mg</sub> halfway between X<sub>Mg</sub><sup>left</sup> and X<sub>Mg</sub><sup>right</sup>. Implicit in the above equations is the assumption that diffusion occurs in an infinitely long solid. This assumption is valid as long as the profiles are long compared with the length affected by diffusion [approximated by √Dt).

To pick promising crystals for diffusion modelling we first acquired low-resolution BSE images with short acquisition times (~60–90 s). High-resolution images were later made across the zone boundaries, with acquisition times of ~30 min. High-resolution images of selected weakly zoned crystals were also taken, so that any minor compositional changes could also be modelled. Profiles perpendicular to zone boundaries were then plotted and
calibrated using the ImageJ software (Rasband, 1997). The profiles were made up of averages across a strip 20–50 pixels wide. The profiles were then fitted by eye in a spreadsheet (LibreOffice Calc), allowing *t*, the diffusion time, to be extracted from equation (19).

## Uncertainties

Temperature has an effect on the opx residence times that can be calculated from the activation energy of Fe–Mg diffusion and the analytical solution to the diffusion equation [equation (19)]. The relationship between the time *t* and temperature *T* is given by

$$t \propto \frac{1}{D} \propto \exp \left( \frac{240000}{RT} \right). \qquad (20)$$

The residence times quoted here are calculated at 880°C; using a temperature of 855°C and the diffusion coefficient of Allan et al. (2013) increases the calculated residence times by a factor of 1.4, whereas calculations done at 905°C decrease the residence times by a factor of 0.7. Using the diffusion coefficient of Dohmen et al. (2016) at 855 and 905°C increases or decreases the residence times by a factor of two.

The oxygen fugacity also has an effect on the diffusion coefficient. The relationship between the times and $fO_2$ is given by

$$t \propto \frac{1}{D} \propto (fO_2)^{-1/6}. \qquad (21)$$

Pre-eruptive oxygen fugacities calculated from Fe–Ti oxides in the Cape Riva and two of the three analysed Therasia pumice fall deposits have a range of ±0.3 log units around the FMQ oxygen buffer. This translates to a difference in near-liquidus residence times of ±12% using the diffusion coefficient of Allan et al. (2013), or ±4% using the diffusion coefficient of Dohmen et al. (2016), relative to the near-liquidus residence times presented here, calculated at the FMQ buffer.

For many crystals, the NLRT estimated for the inner, zone 3-zone 2 boundary are shorter than those estimated for their outer, zone 2-zone 1 boundary (e.g. Fig. 14d and e). This cannot be explained by changes in diffusion coefficient (through changes in *T* or $fO_2$, for example), as any changes that affected the outer boundary would also affect the inner boundary. Resorption also cannot explain this inconsistency, as resorption must have stopped before the growth of zones 2 and 1. The discrepancy could be explained by a mixture of growth and diffusion: if opx growth was more rapid—or the change in magma composition more gradual—during the transition from zone 2 to zone 1 than the transition from zone 3 to zone 2, then the outer boundary would initially be more diffuse than the inner boundary. However, the opposite is seen in the two Al maps of crystals with three zones. The inner boundaries have slightly more diffuse Al than the outer boundaries. Modelling the zone 2 to zone 1 boundary as taking place

in an infinitely long solid may also explain this discrepancy: the crystal rim is close to the zone boundary, and would probably have a composition fixed by exchange with the melt. Using a fixed composition for the rim would reduce the NLRT estimate. However, we cannot rule out the possibility that zone 2 grew first as a skeletal crystal, and that zones 1 and 3 then grew outwards and inwards, respectively, from that initial morphology. A similar growth pattern has been suggested for olivine crystals and plagioclase microlites (Martel & Schmidt, 2003; Welsch et al., 2013). Regardless of the order in which the zones grew, the diffusion profiles provide an upper limit on the NLRT of the opx.

Most of the Fe–Mg zone boundaries can be fitted by a single step function; however, in a few crystals the boundary between zone 2 and zone 1 occurs in two steps (Fig. 14d and e; Supplementary Data Figs S2.13, S2.15 and S2.33). Two steps are often seen in the Al zoning (e.g. Fig. 6c), even in crystals where the zone 2-zone 1 boundary occurs as a single step. This suggests that some of the one-step boundaries between zones 2 and 1 seen in Fe–Mg zoning initially had two steps. The two steps could have merged as they broadened by diffusion. Modelling these boundaries as a single step would give anomalously long NLRT.

Finally, diffusion was modelled at a constant temperature of 880°C. However, zone 2 probably grew from a less evolved melt than zones 1 and 3. This less evolved melt would have been hotter, although the small variation in Mg# (<10 mol %) suggests that the difference would not have been large. If we assume that the diffusion timescale from the boundary between zones 1 and 2 represents the residence time of the whole crystal after the growth of zone 1, we can account for the diffusion that occurred across the inner boundary between zones 3 and 2 during this time. We can then use a different temperature to estimate the residence time of the crystal after the growth of zone 2, but before the growth of zone 1 (Petrone et al., 2016). However, this approach will not work for many of our crystals, as the apparent residence time for the outer boundary is longer than the apparent residence time of the inner boundary (for the reasons discussed above). For those crystals where a two-stage model is possible, using a hotter temperature does not greatly reduce the residence times of the inner boundary. For example, the total residence time estimated for the boundary between zone 3 and 2 in GS11-34a PYX59 profile 1 is 1.5 years at 880°C (Table 7). Using a temperature of 905°C for the first stage (after the growth of zone 2 but before the growth of zone 1) reduces this time to 1.1 years. For GS11-34a PYX63 profile 2 this change is from 0.6 years to 0.4 years. Given that there is uncertainty in how much of the apparent NLRT of the outer boundary is due to growth, and given that we are interested in maximum NLRT, we only use one-stage models for our opx crystals. Because higher temperatures lead to faster diffusion, modelling the crystals at 880°C for the duration of their pre-eruptive storage again provides a maximum estimate of their NLRT.