Earth and Planetary Science Letters Magmatic crystal records in time, space, and process, causatively linked with volcanic unrest

How a has behaved its past is a guide to its future behaviour. Detailed knowledge of what preceded eruptions from speciﬁc volcanoes, and how this can be recognised in real-time, are pivotal questions of this ﬁeld. Here, the physical history of the magma that erupted in 2010 from the ﬂank of Eyjafjallajökull volcano, Iceland, is reconstructed in absolute time and space using only chemical records from erupted crystals. The details of this reconstruction include the number of magma bodies, their geometry, their depth, their relative inﬂation rate and changes to all of the aforementioned through time. Petrology and geodesy (data gathered in real-time) arrive at the same set of conclusions. As such, we report detailed agreement, which demonstrates a causative link between knowledge determined post-eruption via a physical–chemical perspective and knowledge gained syn-eruption from monitoring signals. The composition of olivine crystal cores ( ∼ Fo74–87), and that of the chemical zonation around each core caused by disequilibrium processes, are shown to form systematic patterns at the population scale. Reverse zonation (toward Mg rich) exhibits a constant chemical offset from its crystal core ( ≤ 2 mol % Fo), while normal zonation (toward Fe rich) converges to a single composition ( ∼ Fo75). Conventional petrological models — for instance multiple-magma-mixing across a range of crustal depths — can explain the presence of a range of crystal core composition in the erupted rocks, but cannot explain these patterns of crystal disequilibria. Instead, we describe how a single primitive melt produces crystals over a wide range in composition and generates systematic disequilibrium. Cooling causes crystal production from both roof and ﬂoor of a horizontal magma geometry. Crystal settling causes asymmetric thermal – and therefore compositional – stratiﬁcation of the melt due to progressive insulation via development of a crystal mush at the ﬂoor, a process we term “Crystal Rain”. Crucially, each crystal’s record is both a cause and effect of the internal process of simultaneous fractional crystallisation and settling; no external processes or materials are required. We then extract temporal information from our crystals using Fe–Mg interdiffusion modelling, and combine it with the composition and zonation data. The concept of Crystal Rain is applied, and resolves two thin (metres) sills which are staggered in time and depth, and exhibit different inﬂation rates. Since the approach of integrating crystal chronology within a causative physical framework may be applied to entire volcanic successions, it has potential to yield valuable insights to past, and by inference future, magmatic and volcanic behaviours by deterministic means.


Introduction
Approximately 800 million people live within 100 km of one of Earth's 1551 active volcanoes, ∼80% of which are not adequately Fig. 1. Location map and summary of data gathered in real time. Fimmvörðuháls eruption site (white star) and lateral extent of dilation envelopes of the volcano plumbing system as interpreted from ground deformation and earthquake signals collected in real-time , presented on top of a digital elevation model. The dyke is inclined to the south from the surface to the sills (∼5 km depth; scale bar at the top indicates vertical inflation). Arrows from GPS stations SKOG, THEY and STE1&2 denote the relative magnitude of displacement in the period from November 2009 to onset of eruption on 20th March 2010 in the directions reported ; colours are the same used in later figures. Some figure elements via pers. comm. A. Hooper. are not associated to any lasting record, the deeper past cannot be accessed to aid empirical comparisons and deterministic forecasts (Schmid et al., 2012).
Petrology, by contrast, can read the entirety of the available rock record. Recent research has demonstrated that temporal information using crystal diffusion chronometry (Bouvet de Maisonneuve et al., 2016;Hartley et al., 2016;Kahl et al., 2014;Rae et al., 2016;Ruprecht and Plank, 2013;Viccaro et al., 2016;Zellmer et al., 2005), can be compared on the same time scales as seismicity (Kahl et al., 2011;Longpré et al., 2014), ground deformation , gas sampling (Kahl et al., 2013), and even eruption (Charlier et al., 2012). The important limitations are that such studies can only be conducted after eruption, and crystal chronology is only possible if those rocks contain a tractable mineralogy.
The longstanding challenges have been to 1) accurately determine the dynamics of magma plumbing systems in absolute space and time from volcanic products, and 2) identify direct, causative links between what crystal cargoes record, and any observable signals recorded at the surface. Merging the sub-disciplines first requires us to recognise the different expressions of the same processes within independent datasets, and to link them via shared mechanisms. To achieve this using petrology, each crystal's information must be placed into a single, consistent, context. In order to link these small portions of information into such a context we must conduct investigations at the population scale, for which we require significant datasets. To provide an adequate test, it is necessary to attempt this work where high quality geophysical data are also available.
The 20th March-12th April 2010 flank eruption of Eyjafjallajökull, Iceland, situated in the Fimmvörðuháls pass ( Fig. 1), presents an exceptional case study. With known eruption dates for individual samples Sigmarsson et al., 2011;Thordarson et al., 2011) we can correct timescales of olivine Fe-Mg inter-diffusion (Costa and Morgan, 2010;Kahl et al., 2014) to an absolute chronology, for comparison to pre-and syn-eruption monitoring signals Tarasewicz et al., 2012). Here we first seek to understand the petrogenesis of the magma(s) involved using chemical, textural, and chronological observations, and then compare our results and interpretation with that derived from the geophysical record. Our overriding question is "can petrology provide detailed information of dynamic processes comparable to that determinable by monitoring signals?"

Nomenclature and methods
The term 'macrocryst' is used here to avoid the genetic connotations inherent in terms such as phenocryst (a crystal which grew free in the same magma it is observed within); antecryst (a 'recycled' crystal that grew in a previous magma of essentially the same system); and xenocryst (a crystal with little to do with the petrogenesis of the magma observed; its presence is due to 'accidental' capture; Davidson et al., 2007).
Phase rules are fundamental to petrologic interpretation and a considerable literature and range of practical tools are available to draw upon (see Spera et al., 2016 andMoyen et al., 2016 for recent introductions). The scope of this study limits us to reiterate a single point; phase diagrams are abstracted from the kinetic limitations that exist in nature, which demonstrates near ubiquitous evidence for non-equilibrium processes, as well as incomplete extent of reactions toward equilibrium (Shore and Fowler, 1996). Since dynamic scenarios are discussed throughout, for clarity we define the following.
When the term magma 'Equilibrium' is used, it is implied that the melt is chemically homogeneous, and all formed crystals are of the same composition (per mineral). An uncapitalised 'equilibrium' describes an end-point for changes that may never be realised, but for which there is driving force via equilibration mechanisms such as growth, resorption and diffusion.
Olivine composition is expressed as the molecular percentage of the forsterite (Mg 2 SiO 4 ) endmember (Fo) in the olivine solid solutions, i.e. (Mg 2+ , Fe 2+ , Mn 2+ , Ni 2+ , Ca 2+ ) 2 SiO 4 . Magmatic olivine is essentially a single solid solution between Mg and Fe endmembers; the Mn (tephroite endmember), is usually present as a minor component and Ni and Ca are common trace cations (Deer et al., 1982). An olivine's Fo content is linked to the Mg# of the surrounding melt (Roeder and Emslie, 1970) where Mg# = Mg/(Mg + Mn + Fe) on an atomic basis. Olivine crystallisation and fractionation preferentially removes Mg, relative to Fe 2+ , from a melt (Roeder and Emslie, 1970), driving Mg# down. Subsequent olivine is thus increasingly Mg-poor, forms at lower temperatures, and is described as more evolved. See also supplementary sections 4 and 5 for notes on the chemical system here, including incorporation of ƒO 2 constraints (Gunnlaugsson et al., 2013). Olivine macrocrysts (>1 mm, n = 236) across four samples from the flank eruption of Eyjafjallajökull (Fimm7, -5, -3, -2 in order of eruption) were mapped with high-resolution Z -contrast backscattered electron (BSE) imaging and electron backscatter diffraction (EBSD) techniques. Core-rim profiles were analysed using wavelengthdispersive electron probe microanalysis (WD-EPMA; see supplementary section 2 and also Pankhurst et al., 2017), and a subset (n = 12) of crystal margins were mapped for Ni and P using WD-EPMA, and Mg, Fe and Si using energy-dispersive spectroscopy (EDS). BSE imaging produces maps of relative electron density; for olivine, brighter pixels indicate a higher concentration of Fe (darker = Mg rich). We calibrated the Fo content of an average pixel value located at a given analysis spot to within ∼1 Fo unit (see supplementary section 2). Using per-image calibrations, profiles up to 100 pixels wide were extracted, which allowed us to capture major element chemical distribution with high spatial and compositional resolution, as well as textural and crystallographic context.
The textures and chemical zonation from each crystal were integrated and compared. We first establish the most likely dynamic petrological scenario that fits these data, which is used to define boundary conditions for the retrieval of diffusion timescales (from Integrated olivine composition and textural data. (a) Core compositions (x-axis) from four samples spanning the eruption are presented using a kernel density estimate which utilises a Gaussian distribution function (Ellison, 2002) and a consistent h value of 1.0 (reflecting the 2σ error we calculate for a single EPMA spot). (b) Raw data from (a) plotted according to zonation direction from core. (c) Colour drape over Backscattered Electron image of typical olivine macrocryst exhibiting reverse, then normal zonation towards the high-Fe rim. Note the importance of extracting a profile perpendicular to the c-axis orientation (indicated with arrow direction of plunge). Diffusion is ∼6× faster in this orientation . Choosing the narrowest profile present also minimises the influence of sectioning effects that can cause the apparent profile distance to be longer than the real length. 103 crystals; see supplementary section 6). We then compare the petrological reconstruction with that of remote sensing data.

Results
Macrocrysts comprise ∼30% of the samples' volume, of which 80-85% are olivine and plagioclase in similar proportions, and 15-20% being clinopyroxene. All three mineral phases occasionally occur as glomerocrysts. The erupted melt (∼60-70%) exhibits a homogeneous Mg# value of ∼0.23 , close to equilibrium with groundmass olivine at Fo71 (Toplis, 2005). In contrast, the range of macrocryst core composition extends over 20 Fo mol%, and all but one of the cores are more forsteritic than groundmass olivine (Fig. 2a, see supplementary Appendix A).
The Fo variation of cores is present within each of the four samples (Fig. 2a). Three Fo population peaks are recognised. These are not narrow peaks exhibited by picrites and some other basalts (e.g. Thomson and Maclennan, 2013), but comparatively broad-based populations spanning >5 Fo units.
The main peak (Fo79) and the most magnesian population (Fo86) are present throughout the eruption, the least magnesian crystals (Fo74) are mainly erupted during the first few days. Crystals from all peaks are euhedral, although some Peak Fo86 crystals show some strong resorbtion, see supplementary Fig. 1. Chemical maps reveal that all crystals exhibit homogeneous cores, and all possess chemically-zoned margins ∼50-200 μm wide, which vary in complexity (Fig. 2b, c, see supplementary Appendix B).
An olivine crystal whose composition trends from Mg to Fe-rich toward the rim is said to be "normal zoned". The converse is "reverse zoned". All the crystals except one (whose core is much more evolved than the erupted melt) exhibit overall normal zonation in their outermost margins (usually ∼50 to 150 μm) and <5 μm wide high-Fe rim. Roughly half the crystals of core Fo∼79 (Fig. 2a, b) and ∼85% of those of core Fo∼74, however, exhibit a reverse zone immediately outboard of the crystal core and inboard of the normal zoned margin (Fig. 2).
Qualitative phosphorus and nickel element maps ( Fig. 3) complement the Fe and Mg profiles and BSE images. Of particular note is the consistent co-location of high P zones immediately adjacent to the rimward side of the reverse zone peaks (see also supplementary Fig. 2). Any relationship between high P zones and the normal zones is more difficult to observe. Ni corresponds spatially with Mg in both reverse and normal zoned crystals.
When we compare the core composition of olivine within reverse and normal zones of the same crystals, two different arrays are observed (Fig. 4, see also supplementary Appendix A). In the majority of cases the compositions of reverse zones are ≤2.5 Fo units more forsteritic than the core. A strong positive array is formed on graphs of Core XFo and Reverse Zone XFo, which spans ∼10 Fo units (Fig. 4a). Normal zone profile shapes include a strong linear component (supplementary Fig. 3), taken to indicate crystal growth rather than purely diffusion, which would lead to sigmoidal or continuous curves depending on the boundary conditions. Some overall curvature implies a sub-dominant diffusion component, or simply evolving melt conditions. By contrast, reverse zones are best described as displaying classical, asymmetric diffusion profile shapes (Fig. 2c).
Modelling Fe-Mg profiles as pure diffusion in the region of the reverse zones returns very good fits (supplementary Fig. 4), although a small number did not pass stringent quality assurance/quality control (QA/QC) and were discarded (see supplementary section 6 for details). Calculated diffusion timescales (using the above boundary conditions, and the same numerical methods as previous work; Hartley et al., 2016;Rae et al., 2016) range from a few hours to a few months. There are proportionally more crystals with shorter timescales than longer timescales, which we illustrate on a cumulative frequency diagram (Fig. 5).
There is no obvious clustering, but rather a progressively increasing number of reverse zone initiation events through time until late February, when this trend flattens (Fig. 5). In early March, the cumulative frequency suddenly steepens again, and after the eruption commences it tails off (which is an artefact of fewer In both cases (normal, a; reverse, b), the c-axis is at a high angle to the sectioned plane. Black bars (in aii and bii) are 50 μm. Each map is a linear colour drape scaled across image intensity range. Note the progressive growth in the normal zoned crystal (a) that matches Ni, Fe, Mg patterns. By contrast, in the reverse zoned crystal (b), Ni, Fe and Mg have diffusively relaxed from a P-defined crystal edge (T2 in bv) into the crystal interior, best illustrated by the rounded interior of the reverse zone with respect to the pointed exterior. Time series (black, green, blue, red) are relative to each crystal; no equivalence to each other or other time-series in this study is intended, and are presented to illustrate that wherever non-linear chemical gradients exist, diffusive relaxation is expected, which can be simultaneous with further growth outboard. See also supplementary Fig. 5. The colour overlay across each 8-bit image was chosen to highlight chemical contrast in each case; no equivalence of concentration between maps is implied. samples contributing to the frequency distribution as the eruption progresses past the eruption dates of analysed samples). Individual samples contain crystals which produce the same curve when plotted together or apart.

Discussion
Our olivine chemical population data (n = 236) are consistent with Moune et al. (2012) and Viccaro et al. (2016) who used 28 crystals each. A simple melt evolution can be approximated involving mantle-derived magma cooling down an olivine liquidus, joined later by plagioclase and clinopyroxene (see supplementary sections 4 and 5). Our timescale data are consistent with those of Viccaro et al. (2016) who used Fe-Mg diffusion profiles from 9 olivine crystals from a lava flow produced late in the flank eruption. With our data density, fully integrated composition, texture, zonation directions and diffusion timescales (>10× that of Viccaro et al., 2016), we can attempt to describe where and when each crystal originated, what happened to each after their formation, and their relationships with the plumbing system.

Olivine petrogenesis and relationship to physical processes
Some crystals record high frequency Fe-Mg oscillations in an overall normal zoned margin occupying the final few 10 s of μm from the rim. These must be formed in the final few hours before quench in order to be preserved (and not diffusively relaxed), which may contain insights to eruption dynamics. The potential for other crystals to have formed this texture, and have relaxed at least to some degree precludes detailed comparison here. Furthermore, they represent far shorter timescales than those associated with the geodetic signals of unrest.
By focusing on the compositions of the homogeneous crystal cores, and those of the adjacent zones, we consider each crystal from its local Equilibrium origin forwards in time. A reverse zone immediately outside the core occurs in about half the crystals and represents a traceable, comparable point in those crystals' history. A logical explanation has to cover why reverse zonation is present in only ∼half of the crystals, as well as explaining how a systematic compositional offset (of <2.5 mol%) was formed (hereafter 'reverse offset').
First, we used trace element distributions to check the relationship between Fe-Mg and textural and chemical histories. Assuming P in melt was not saturated, the inner edge of high P zones within olivine picks out crystal edge boundaries prior to bursts of growth. Hence, it can be used to chart crystal texture with time (Milman- Barris et al., 2008). This is because P has both low com- A continuous array from several months before eruption is observed, which continues after the onset of eruption. Cumulative frequency is presented by plotting the start date on the x-axis (calculated simply by taking the sample eruption date and subtracting the diffusion timescale) and spreading the data equally in the y axis according to rank order of appearance forward through time. The upward concave shape reflects a weighting in the data to shorter timescales. A brief hiatus of start dates occurs between late February and early March. It divides an older, broad, concave-upwards curve (<5% of the crystal timescales, including their uncertainties, are older than the time window chosen for display) from a younger, steeper curve that continues to shortly after eruption onset. Horizontal lines show confidence limits, per crystal, equivalent to a 0.2 log 10 unit uncertainty, reflecting a temperature uncertainty of ±15 • C. A difference of 0.5 log units on fO 2 would correspond to a time difference of 21% absolute, yet the available Fe speciation measurements (Gunnlaugsson et al., 2013) suggests that the fO 2 difference along the olivine liquidus that produced the macrocrysts was less than this (see supplementary sections 5-6). Note that the calculations have smaller uncertainties at shorter timescales; uncertainties are also asymmetric when displayed on a linear time axis. The uncertainty envelope is spiky as uncertainties vary in proportion to the time of sample quench (eruption), meaning that the exact time of eruption of a given sample (known for our sample set) is important in constraining the uncertainty for crystals derived from that sample. patibility and low diffusivity in the olivine structure (e.g. Watson et al., 2015). In contrast, Ni is highly compatible, and its diffusivity is comparable to that of the Fe-Mg binary interdiffusion. This implies that Ni will co-vary with Mg during growth and diffusively equilibrate within crystals on similar timescales to Fe-Mg (Costa and Morgan, 2010). P maps reveals a series of discrete growth stages within all crystals ( Fig. 3 and supplementary Fig. 2) as indicated by the presence of high-P zones (Milman-Barris et al., 2008). Ni maps mirror those of the Mg distribution. Normal zones are most simply explained by progressive growth during melt evolution (to higher Fe and lower Mg) after crystal nucleation in Equilibrium conditions (homogeneous core).
Of particular note is the distribution of P and Ni in reverse zones. In each case there is a clear high-P boundary -an earlier crystal edge -which is now located slightly rimward (by <10 μm) of the Mg (and Ni) peaks (Fig. 3). The simplest explanation is that each of these crystals' Equilibrium growth was stopped, melt conditions changed to a higher Ni concentration for a period of time (reverse zone), before growth recommenced along a fractional crystallisation trend (normal zone). Fe-Mg interdiffusion has seemingly taken place inside the crystal with invariant boundary conditions for much of the duration (Fig. 2c). This implies crystal cores formed in one regime, before residing in another that was slightly hotter, yet physically and chemically stable. We now discuss three types of magmatic scenario, and consider if each may bring about this crystal-melt behaviour.

Scenario #1, magma heating
The simplest physical scenario (Scenario #1) is one of heating up an Equilibrated, pre-existing magma; phase diagrams predict a shift to higher Fo equilibrium. Reverse zones in Fo could be ex-pected to form, but dissolution is also predicted. In addition to simplicity, scenario 1 appeals because it requires only one melt composition, which is what we observe.
Scenario #1 does, however, suffer several difficulties against our observations. The reverse zones are straight-sided (except on corners; Fig. 3), and lack evidence for dissolution. Resorption margins are only observed in the Fo86 peak (the most primitive), which suggests these features formed by those crystals responding to the cooler, and more evolved conditions of the final carrier melt.
An extension of Scenario #1 is injection of hot, primitive melt into a pre-existing cooler and more evolved magma body. This is a prima facie satisfying idea, especially so for olivine, because it provides both heat and higher Mg# conditions to the system, as well as replenishing Ni. However, this does not explain the wide range of Fo in crystal cores, for on a heating trajectory, olivine loss, not formation, should occur. It also is inconsistent with our observation of a uniform and evolved melt composition at eruption.

Scenario 2, magma mixing and hybridisation
Should two different magmas pre-exist, each containing their own crystals in Equilibrium (Scenario #2) we could describe a simple process of mixing. Since the new equilibrium would lie somewhere between the two, each magmas' crystals will re-equilibrate toward a hybridised 'mid-point'. Catching the process in the act would mean observing two directions of zonation linked to crystal core composition; primitive crystals will become normal zoned, and evolved crystals will be reverse zoned (and/or undergo resorption).
Scenario #2 does not predict crystals at the same core composition, of which some are normal zoned, and others reverse. Our data contains such evidence in abundance. This complication might intuitively be overcome by appealing to a more complicated mixing history, whereby subsequent pulses of the fresh, primitive magma 'catch' the system already in a mixed state, but prior to full crystal equilibration. Yet this describes mixing along the same liquidus as defined by the initial mixing, and does not resolve the bi-directional zoning from a single core composition problem. It also reduces interpretation to single-crystals; since each have nonunique solutions, this approach cannot support meaning attached to any patterns observed at the crystal population scale.
Perhaps the most compelling observation that argues against batch magma-mixing is the consistent size of offset between reverse zones and the cores. Melt of highly specific composition only comes into contact with highly specific composition olivine cores, right across the range. Since the simplest magma mixing scenario does not explain this observation, and increasing complexity only serves to decrease the likelihood of such a systematic pattern, any form of hybridisation involving distinct magma bodies must be so highly regulated that it would seem implausible.

Scenario 3: incomplete cross-mixing and parallel evolution
Crystals from the Fo74 and Fo79 peaks must be contacted by melts both more and less evolved. Scenario #3 is to explicitly allow the introduction or production of crystals across a range of core compositions in parallel, by coincident magmas and/or involvement of host rocks. This would allow 'cross' mixing (crystal transfer) to occur, where crystals of the same composition come to be surrounded by either evolved (perhaps settled into a pocket of evolved magma?) or primitive (a new injection?) melt. While the concept is at best stochastic, it forces us to consider that the crystals and melts can become decoupled and behave semiindependently.
Nature indeed contains abundant evidence for complex, opensystem scenarios (Martin et al., 2010;McLeod et al., 2012;Shore and Fowler, 1996;Streck, 2008;Tepley III et al., 1999), across a range of magma compositions (Bouvet de Maisonneuve et al., 2015;Dobosi and Fodor, 1992;Pankhurst et al., 2011;Vernon, 1990) and timescales (Bindeman and Melnik, 2016;Martin et al., 2008;Morgan et al., 2004;Rioux et al., 2016). While often the most reasonable to invoke, Scenario #3 is the one most removed from the predictive qualities of phase diagrams. Thus, as complex explanation can fit the complexities of natural crystal data, Scenario #3-like models have little in the way predictive capacity because there are multiple plausible solutions for each crystal and for the overall petrogenesis.
Specifically, Scenario #3 does not provide for highly systematic trends and patterns at the crystal population level, whereas we observe two systematic trends composed of all but one crystal (Fig. 4). While crystal transfer must occur (since Equilibrium, or fractional-crystallisation alone or acting together cannot explain our patterns), it does not appear to be stochastic.

Toward an explanation for the Fimmvörðuháls magma
Any physical magmatic model must explain both the melt homogeneity, and the heterogeneous crystal core chemistry. It must also explain the limited range in reverse zonation magnitude with respect to the crystal core, yet allow a larger range in normal zonation magnitude from cores to be produced, and a convergence towards a narrow composition range.
Fo∼86 crystals are in equilibrium with mantle-derived magma (Viccaro et al., 2016; see also Thy, 1983 andKress andCarmichael, 1991), as opposed to being xenocrystic mantle olivine. This implies that mantle-derived magma was involved. The erupted melt composition was far more evolved than that in Equilibrium with Fo86. Together, these observations could be explained by a mantlederived magma cooling and crystallising in the crust, without needing a second magma to mix with, which is consistent with melt inclusion data Moune et al., 2012; see supplementary sections 4 and 5).
This resembles none of the three scenarios above, and more a fractional crystallisation scenario where cooling leads to internal differentiation (Couch et al., 2001) via the production of crystals along a liquidus. It does not address the presence of reverse zones, however, or why roughly half of the crystal population exhibit them, or why there are three distinct population peaks in core composition.
To a first order, crystals nucleate when and where a melt is undercooled (Roeder and Emslie, 1970). Zonation is subsequently imposed, indicating either a change in the local environment, or that the crystal has moved from where it first grew. Absence of zoning implies nothing has caused a departure from an Equilibrium state. Working backwards from the rim, at their latest stage, all crystals appear to have been bathed in the same melt in equilibrium with ∼Fo75.5. Using a kinematic approach (Kahl et al., 2014), this represents their penultimate environment en route to eruption. Using the broader chemical system, the most likely depth for this environment is ∼4.5-5.0 km (Viccaro et al., 2016). This leaves the reverse zoned crystals to explain, and why there are population peaks.

A crystal rain; dynamic equilibrium in a spatial and temporal, in addition to PTX, context
Macrocrysts with reverse zones have core compositions from ∼Fo73-83 (Fig. 4a), and are within the range of core compositions of only normally zoned crystals (Fig. 2b). The working hypothesis involving a single magma limits the range of conditions that a single olivine composition can crystallise. We must consider that the reverse zonation is due to each crystal departing their initial, lo-cal environment, and migrating into one slightly, but consistently, more primitive. Why, and how?
Both questions can be resolved by considering the divergent histories a crystal may experience after originating either in the upper or lower thermal boundary layer of a cooling sill (Marsh, 2015), see Fig. 6. Olivine is ∼25% denser than an equilibrium melt containing ∼1.0 wt% water (Ochs and Lange, 1999). The low calculated melt viscosity of 20-40 Pas (Giordano et al., 2008), together with the dimensions of crystals (all >1 mm) would allow efficient crystal settling (Holness et al., 2017;Martin and Nokes, 1988).
Sedimentation changes the dynamic thermal structure of a sill (Marsh, 2015). A mush pile of olivine crystals accumulating at the floor (Bergantz and Ni, 1999) causes insulation of the lower cooling surface (Marsh, 2015). In contrast, the crystal settling (Bergantz and Ni, 1999;Martin and Nokes, 1988) maintains narrow isotherms at the roof, and thus a dynamic, asymmetric thermal profile is set up whereby greater undercooling is experienced at the roof of the sill with respect to at the floor.
The range of melt composition and temperature present at each moment through time can be visualised as an arc of the olivine liquidus (Roeder and Emslie, 1970), rather than a single point (supplementary Fig. 5). This arc produces crystals of different compositions, across the vertical thickness of the sill through time, focused in the thermal boundary layers (Fig. 6). As the sill loses heat and crystallises, the arc migrates along the liquidus, promoting the generation of a population peak (supplementary Fig. 5).
Roof-zone crystals are produced at one position on the liquidus, and by departing the cooler thermal boundary layer and raining onto the mush pile, they become bathed in more primitive and warmer melt; a higher position on the liquidus (i.e. Kohn et al., 1989). Thus they cease growth, and begin to diffusively reequilibrate, forming the reverse zone and its characteristic profile shape (essentially with boundary conditions of an external buffer and no overgrowth, as compared to a square wave function within a crystal due to an overgrowth event). The maximum magnitude of the reverse zonation is determined by the maximum temperature contrast across the sill at any point in time (T1-4 in Fig. 6). As such, there is never an opportunity in time and space for highlyevolved crystals to come into contact with highly-primitive melt.
Crystals formed directly on, or close to, the floor do not move; and so no reverse zone forms. At any one point in time, the overlying melt is always more evolved than all olivine crystals in the mush except the most recently formed crystals derived from the roof zone. In this mush, crystals are volumetrically dominant and their packing limits the potential for chemical diffusion within the interstitial melt. When the mush is mobilised before eruption, all crystals become bathed in the residual cooler, more Fe-rich and volumetrically plentiful melt; what little trapped melt was present in the mush is quickly mixed and normal zones are formed on all crystals.
Thus a chemically-stratified mush comprised of both equilibrium, solely normal zoned (floor-derived) and reverse then-normal zoned (roof-derived) crystals (Fig. 6) explains how cores of the same composition can have divergent zonation senses (Fig. 2b); they reflect divergent formation histories. Of critical note is that a wide range of crystal-melt pairs form in an extremely limited vertical extent, as an inherent part of the systems' evolution towards equilibrium.
Importantly, there is no whole scale cross-mixing required, and the many endmembers that the batch magma paradigm contains are avoided. This is a continuum model, where crystal cores are produced in an ordered fashion, consistent with a single liquid line of descent. The production and migration of these crystals themselves define the processes that generate the subsequent zonation patterns. This is in contrast to the intuitive philosophy that a continuum process should produce simple zoning (and explicitly not Fig. 6. Crystal Rain model. Crystal Rain is a shorthand description of the process by which differential nucleation and growth inside a magma body, caused by crystal settling that imposes an asymmetric thermal and compositional gradient, incorporates crystals into an insulating mush horizon according to dynamic, yet systematic chemical stratification in terms of both crystals and local melt. Subsequent, possibly repeated, entrainment into overlying magma could produce oscillating normal zones around all crystals. account for reversals or discontinuities; Bouvet de Maisonneuve et al., 2016).
Crystal Rain is an effort to rationalise both Equilibrium and fractional crystallisation where neither can alone explain our data, and includes short-distance crystal transfer. Part of the novel distinction is that in order to explain the patterns, all three classic mechanisms must operate in a reasonably self-regulated manner. To develop the Crystal Rain scenario into a testable model with predictive capacity, extensive thermochemical and physicochemical work needs to be completed. What can be achieved here, however, is testing whether the petrologic scenario is consistent with chronological data. At a crystal scale, these data are independent from the way phase rules underpin the classic mechanisms, since phase rules are atemporal. At a magma plumbing system scale, we can turn to remote sensing records.

Diffusion timescales in the context of Crystal Rain and mobilisation
Local boundary conditions remain relatively stable within a mush, certainly over the comparatively short timescales we observe (e.g. Thomson and Maclennan, 2013). The simple reverse zone profile shapes and the required boundary conditions to model them point to a stable physical environment for the majority of the crystal's history. A sill's horizontal geometry lessens the potential for mushy portions of a magma body to experience gravitational collapse after forming.
Crystal diffusion timescales represent the duration spent at high temperature following establishment of disequilibrium . A non-convecting mush retains heat with respect to the overlying convecting melt. Thus the crystal 'clocks' 'tick' in the mush at comparable speeds due to the similar temperature .
The duration between mobilisation of the mush and eruption (quench) is implicitly included in these calculations. The normal zone forms by overgrowth during this time, and these conditions are slightly cooler than that of the mush. Now, the outer boundary conditions have changed; where once was a more primitive melt than the core is now more evolved olivine than the core. This is unlikely to be a major source of error. Any cooling experienced en-route would be reflected both in the final melt fraction (70%) and volcanic behaviour (fire-fountaining), which are both consistent with near-magmatic temperatures on eruption, which means the temperature (diffusion speed) is still a good approximation. In addition, the proportion of time crystals experience en-route is very small compared with that in the mush, which is consistent with interpretations of combined elemental and isotopic profiles in natural olivine crystals (Sio and Dauphas, 2017). If the proportion of time crystals experience en-route were larger, we would expect to see a pronounced decoupling between P and Ni zonation and the Fe-Mg profiles, since the reverse zone "well" would be "filled up" (with Fe). Instead what we observe is a good match between the early rim position (see Fig. 3) and the location of the profile max Fo. Sequential modelling (e.g. Kahl et al., 2011) makes no significant difference to the final timescale.

Plumbing system behaviour in the lead-up to eruption at Fimmvörðuháls
The diffusion timescales approximate the time of sedimentation into the mush. This means the shape of a timescale rank order plot (Fig. 5) is a function of the rate of crystal formation, which is in turn the direct result of the rate of magma volume undergoing cooling. The slope of Fig. 5 is thus a qualitative proxy for the rate of energy loss into the crust; a function of the sill surface area as it intrudes.
Since we can model the liquid line of descent in this chemical system, as others have done in great detail (Viccaro et al., 2016), we can ask at what depths and temperatures (and ƒO 2 ) olivine crystals of certain compositions are in Equilibrium. The population peaks provide the only means of assessing where that equilibrium is. As olivine is produced on a pathway to an equilibrium point, the minimum Fo is a useful 'check-point' and is the best evidence that be used to parameterise a depth calculation (see supplementary section 5).
We may reconstruct the active plumbing system in the lead-up to eruption, as illustrated in Fig. 7. Stage (1) involves the ascent of mantle magma which forms a sill in the crust (2). The two cool- 2) Crystal-poor parental magma ascends to a level of neutral buoyancy and/or zone of stratiform weakness where it begins cooling and crystallising down an olivine liquidus towards an equilibrium at lower pressure and temperature (∼5 km depth: see text and Supplement). 3) Crystal Rain initiates and builds a chemically stratified mush pile. Reverse zones form around roof-derived crystals (T1-4 occurs continuously in this part of the plumbing system; see Fig. 6). 4) Mobilisation into overlying cooler, more evolved, melt causes growth around all crystals. Position of GPS stations added for reference frame (see Fig. 1) and to define lateral extent of the sills (not their depth).
ing surfaces induce Crystal Rain (3) and the Fo79 peak forms. The Fo74 peak can be explained most simply by suggesting the same mechanism that formed the Fo79 peak also formed the Fo74 by initiation of a second sill, probably a little shallower (4).
This follows the reasoning that since cooling drives the crystallisation, the initial thermal environment in the crust must play some role in defining where the equilibrium point would be for a given magma composition. In the absence of more sophisticated, multiphase modelling, we may simply assume the crustal geotherm provides cooler conditions at shallower levels. It follows that the equilibrium the magma evolves to would be further down the olivine liquidus.
With a physical scenario that is self-consistent, we now integrate timescale data with these concepts. In Fig. 8a we show that the two time periods observable in Fig. 5 also coincide remarkably with a drop in minimum Fo content of olivine cores. This is best explained by representing two different periods in time and chemical space (Fig. 8b). We interpret these observations as reflecting the intrusion of two sills that were staggered in time (Fig. 8c).
This temporal staggering has implications for the magnitude of the population peaks. We may now view them as the production of crystals, from ∼Fo83 to ∼Fo76 in Sill #1, and from ∼Fo83 to ∼Fo73 in Sill #2, and thus the Fo79 peak is the tallest because it reflects the super-positioning of one sill's journey to equilibrium upon the other. With a constant influx of primitive magma that rapidly evolves to its equilibrium temperature due to the large cooling surfaces, a 'bank' of crystals is formed across a range of compositions, which subsequently are mobilised and withdrawn from the sills. The increasing primitive peak in olivine core composition likely reflects progressive skimming off the top of the crystal mush as the eruption progressed, working downwards where a greater proportion of primitive olivine resided.

Reconciliation with monitoring
Crystal timescales in period 1 are a match with the period of initial inflation observed by InSAR (Fig. 8d, e; Table 1; Sigmundsson et al., 2010). The match is further, and independently, substantiated in the depth and geometry of each sill. For example, the rapid rise in the slope of the crystal data coincides with the sudden in- Fig. 8. Results, interpretations, and geophysical data for comparison. (a) Crystal diffusion timescales are corrected for time of eruption and plotted in rank order on a true timeline. Crystals are colour coded by composition and symbol-coded for sample. The grey envelope represents uncertainties as in Fig. 5. (b) First-order interpretation of integrated timescales and crystal core compositions reveal two distinct periods, with a short hiatus in crystal production detected between the two. (c) As viewed through the Crystal Rain model, these distinct time periods are most simply interpreted as reflecting the development of different sills; the steepness of the cumulative frequency in (a) could be related to propagation/crystal production rate and the minimum olivine composition corresponds to an equilibrium point, in part reflecting pressure (depth). (d) Selection of ground-based geophysical data; 3 GPS stations SKOG, THEY and STE1/STE2 and the envelope of detected earthquakes each day (grey) from , whose interpretation is summarised in (e), which is a close match with (c).

Table 1
Time-linked observations of magma plumbing system dynamics as determined independently by remote sensing and petrology.

Absolute time Observation and rationale for interpretation
Remote sensing  Petrology ( crease in inflation seen by InSAR, and with increased earthquake rates. The geophysical data alone describe the lateral position of each sill (Fig. 1), whereas the petrological data alone describe the chemical evolution and intensive parameters of pressure and temperature as it evolves within each sill. In each case, these parameters help to reconstruct an equivalent and detailed history of the pre-eruptive behaviour from petrology, as compared to that from monitoring (see Table 1, and compare Fig. 8c with e). Seismicity 10-15 yr earlier, at deeper levels (5-6 km), likely relates to a similar style of magma intrusion . It is unclear whether those thin sills remained connected and could have contributed to the erupted material, or cooled to some temperature (or indeed completely froze) to become mechanically indistinct from country rock. Since a) the year 2010 sills developed at shallower levels than those at least a decade earlier and b) a self-consistent, 3D plus time petrological model can be described without requiring any pre-existing magma (Fig. 7), it is simpler to view the earlier seismicity and intrusions as relating to the underlying mantle source beginning to release magma, but cooled in the crust prior to the 2010 pulse.

Conclusions
It is possible to explain the variety of compositions and textures in olivine from the 2010 eruption at Fimmvörðuháls by a continuum model of magma migration into the crust, where relatively short-term crystal fractionation, sinking and storage occurs in sills. Crystal Rain affords details of the petrogenesis to be recovered such that reconstruction of the magmas dynamic physical and chemical behaviour within these sills, and in absolute time and space can be achieved. This study demonstrates that petrologic, geophysical and remote sensing interpretations of preeruptive plumbing system behaviour can independently arrive at the same causative explanations (see Table 1). We conclude that it is plausible to produce a specific and detailed volcano monitoring record using petrology; a causative proxy for the kind of information gained in real-time.

Petrologic pressure estimates
Geobarometry uses knowledge of pressure-sensitive elemental partition coefficients to retrieve depth estimates for crystal-liquid equilibrium (Putirka, 2008). Crystal Rain describes a dynamic system that is not in equilibrium, but evolving towards it, produc-ing and containing diverse crystal chemical populations simultaneously within a highly-restricted vertical extent. The concept may provide new context to geothermobarometric and melt inclusion studies.

Limiting degrees of freedom in interpretations of complex crystal records
Of special importance are the inherently non-unique solutions generated by scenarios within the batch-magma mixing paradigm. Rather than complicating matters, we show that by integrating and rationalising chemical and temporal data into a dynamic physical scenario, Crystal Rain can in fact simplify the number of possible solutions. In our solution, each crystal is contextualised when viewed as both a cause and effect of a simple process, and this leads to a much simpler inferred system in terms of geometry and the number of magmas involved. Here this concept has allowed us to retrace magma pathway(s), describe storage geometry and depth of cooling and evolving magma bodies, and understanding the timescales involved of each processes using only erupted products.
Future work aims to combine kinetic and kinematic processes during evolution of a sill of different thicknesses, and compare these results with detailed observations and interpretations derived from exposed fossil sills (e.g. Holness et al., 2017 and references therein). Such comparisons, however, are inherently limited since our observations capture processes operating on timescales such that subsequent relaxation of Fe-Mg gradients in olivine are likely to have been erased during any extended, hot crustal residence.
Liquid lines of descent as computed from petrogenetic programs such as MELTS (Ghiorso and Sack, 1995), Petrolog3 (Danyushevsky and Plechov, 2011), and EC-RAFC (Spera and Bohrson, 2002), provide a 'snapshot' Equilibrium view of a system at various times through magmatic processes, but as yet do not accommodate disequilibrium processes. We would encourage the next generation of numerical models to link thermal gradients within cooling magma bodies with the chemistry of saturated solid phases, for this would allow the prediction of dynamic, albeit often subtle, chemical gradients that appear to be detectable in natural data.

Volcano monitoring
Catastrophe models are generated using probabilistic techniques to quantify volcanic risk. These draw upon the distribu-tion, volume and type of extrusive product with the volcanic rock record -the 'back-catalogue' of eruptions -to define the hazards which help gauge likely impacts (Schmidt et al., 2011;Todesco et al., 2002Todesco et al., , 2006. What is required to radically reduce risk, however, are deterministic models of volcanic behaviour with specific predictive capacity (Schmid et al., 2012).
Adding pre-eruptive behaviour to this 'back-catalogue' would give volcano observers opportunities to contextualise signals of unrest, and form the basis for empirical comparison in real time. Reconstructing dynamic magma behaviour in time and space for prior to eruptions across a volcanic systems' lifetime may reveal higher order patterns of what a typical and non-typical preeruption period 'looks like'. Such insights would be of great value to decision makers in hazard mitigation and civil defence. For instance, real-time observations that are not consistent with any previous pre-eruption period could be flagged as either low likelihood of eruption, or unlikely to produce an eruption similar to those previously analysed. While the 2010 Fimmvörðuháls eruption took place while real-time observations were made, these are not necessary to vastly improve forecasting and crisis management. These data cannot yet be used to determine the duration of the eruption, yet it is perhaps more important to know what a typical plumbing system behaviour (i.e. inflation at a given rate), its potential eruptive outcome, and its timescale is.
At present there is an epistemic gap between state of the art knowledge in physical volcanology and volcano monitoring required to produce deterministic models. This is due to the great range of time volcanoes are constructed over, and the comparatively small window of time active volcanism has been monitored. This study demonstrates the gap can be overcome, without needing to experience numerous eruptions to link pre-eruptive unrest with the surface outcome (by which time it could arguably be too late).
While proxy volcano monitoring records will inevitably have limitations of timescale resolution, the approach holds potential to allow us to extend monitoring-equivalent data back into the geological record. It is likely that hazardous volcanoes of suitable mineralogy will be the highest priority for this type of assessment, yet this approach can be equally applied to volcanic systems where no current monitoring network is present.

Author contributions
Sample preparation, data acquisition, conception of the petrogenetic scenario and manuscript formation was conducted by MJP under the remit of NERC grant NE/J024554/1 to DJM, who wrote the diffusion modelling algorithm and performed the initial modelling; TT who supplied material, and SL. All authors contributed to discussion of the concepts and writing of the manuscript.

Author information
We declare no competing financial interests. Correspondence and requests for materials should be addressed to mpankhurst @ iter.es.

Data statement
All data will be made available via the NERC (UK) National Geoscience Data Centre, and/or on request from the authors.