The time elapsed between magma recharge and eruption is recorded by the zoning of Group 1 (Fo-like) elements at the rim, whereas the time elapsed between skeletal growth (i.e., just after olivine nucleation) and eruption is recorded by the zoning of Group 3 (hybrid) elements in the core. Both timescales can be quantified by kinetic modeling using a 2D diffusion model applicable to complex geometries. The lattice Boltzmann model described in the “Methods” section has been applied to coupled Fe–Mg zoning near rims and Ti and Sc zoning in cores.

Timescales between mixing with magma recharge and eruption

Timescales between the interaction of mafic recharge with resident magma and eruption were constrained for nine olivine phenocrysts using Mg X-ray maps converted to quantitative Fo maps. Initial conditions consist of a crystal of the shape of the analyzed olivine with a uniform composition corresponding to the plateau measured in the core. Boundary conditions were maintained constant at the measured composition of the rim. Any discrepancies between the measured and the modeled concentrations were computed at each time step by the least square error in concentration integrated over the grain. The Fourier number, Fou, corresponding to the best fit (i.e., minimum least square error) was converted to a timescale using appropriate values for diffusivity. The diffusivities of Fe and Mg in olivine are anisotropic, strongly dependent on temperature and Fo content, and slightly dependent on pressure and oxygen fugacity (Table 3; Dohmen et al. 2007). Pressure and oxygen fugacity were considered to be 0.1 MPa (atmospheric) and 5.10−4 MPa (close to the Ni-NiO buffer at 1130 °C). The temperature at which diffusion occurred is the most difficult parameter to constrain, as (1) it was probably not constant during diffusion (unlike in our model) and (2) different crystals probably experienced different temperature histories depending on whether they were residing in evolved mush or carried by hotter recharge magma. Timescales are therefore reported as a function of temperature.

Table 3 Olivine diffusivities of Fe–Mg, Ti, and Sc used in this study Full size table

Hybrid magma temperatures calculated by olivine–melt equilibria at the crystal rim cluster around 1130–1150 °C in all three units (anhydrous estimates; Bouvet de Maisonneuve et al. 2012). If the average water content in the resident mush was 2 wt% (in agreement with melt inclusion data), pre-eruptive temperatures could have been closer to 1060–1080 °C rather than 1130–1150 °C, after applying the correction of Médard and Grove (2008). If a single, characteristic temperature is to be used to reproduce diffusion under non-isothermal conditions, it should be close to the maximum temperature at which diffusion occurred as diffusion rates vary exponentially with temperature (Chakraborty and Ganguly 1991). We therefore considered diffusion to have occurred at the elevated temperatures of the hybrid magma rather than the cooler temperatures of the mush. Below we report timescale estimates for both anhydrous and hydrous (2 wt% H 2 O) conditions, as these represent two end-member scenarios that bound the range to be expected.

Two-dimensional diffusion modeling reproduces the Fo content of the maps within 1–1.5 mol% (Fig. 7). Larger discrepancies occur at broken edges (e.g., 1957-A3, 1957-B6, LF3-A27), in crystals that show the effect of diffusion in the third dimension (e.g., UF3-D14 and E3), or where melt channels are connected to the exterior of the crystal (e.g., LF3-A27). These effects are small and do not substantially influence the final result, as zoning toward the rims in particular is well reproduced. Model results overlap for the three eruptions, with timescales of a few months to a few years for a temperature range between 1150 and 1060 °C. In more detail, elapsed times between magma recharge and eruption seem to have been shorter prior to the 1957 summit eruption than prior to the early Fissural 3 eruption (LF3, ~1850 AD) at a given temperature. If hydrous pre-eruptive temperatures of 1060–1080 °C are considered, magma recharge took place between 8 months and 3 years prior to the 1957 eruption, between 3 and 4.5 years prior to the UF3 eruption, and between 3 and 9 years prior to the LF3 eruption. With anhydrous pre-eruptive temperatures of 1130–1150 °C, timescales would be reduced by approximately a factor of three.

Fig. 7 2D lattice Boltzmann modeling (LBM) of the Fo zoning in the nine olivine crystals studied from the 1957 (a, d), UF3 (b, e), and LF3 (c, f) tephra deposits. (a-c) EPMA Mg X-ray maps converted to quantitative Fo maps, best fit Fo zoning map from LBM, and spatial distribution of the error (here shown as the difference in concentration between modeled and measured Fo) are shown for all crystals. The color bar on the left of the maps shows the Fo content in the crystal, while the color bar at the base of (a) shows the magnitude of the error for all crystals. d-f Time elapsed between magma recharge and eruption as a function of temperature. The most likely range in temperatures and their associated timescales are shown with gray shaded areas (hydrous) and gray dashed lines (anhydrous estimates) Full size image

Shorter time lags between the 1957 eruption and preceding recharge events are consistent with the long-term eruptive behavior of the main summit cone vent, from which a significant eruption occurred only 10 years prior to the very large 1957 event (Naranjo and Moreno 2005). Eruptions from the strikingly steep-sided main cone vent system, which is controlled by N–S-trending normal faults, are far more frequent and more voluminous on average than those from lower flank vents. The so-called Fissural 3 event (LF3-UF3 tephra) of Naranjo and Moreno (2005), to which we ascribe an age of ~1850 AD, is a large lower flank vent eruption, but the two preceding phases of activity from this vent area are prehistoric and may have occurred up to several 1000 years B.P. We infer that the higher eruptive magma flux from the main summit vent system is associated with a much higher rate of background recharge, hence, a much higher probability that olivine crystals from large summit eruptions will record the higher pre-eruptive temperatures and shorter recharge-eruption lag times.

Timescales between early crystal growth and eruption

Timescales between olivine skeletal growth (i.e., immediately after nucleation) and eruption were constrained for two crystals (1957-A3 and B1) in which trace element zoning in the core was not overlapping with zoning at the rim (i.e., overprinted by subsequent recharge events). Initial conditions consist of a crystal with the shape of the analyzed olivine and a skeletal core of the shape of the high-P zone (Fig. 8). The skeletal core has initial enrichments in Ti and Sc related to the enrichments in P and dictated by the growth rate, partition coefficient, and diffusivity in the melt. The outer part of the crystal has no enrichments in Ti and Sc (i.e., value of 1) and defines the boundary condition. The P-rich oscillatory zoning near the crystal rim was not considered for modeling as (1) we believe that it belongs to a second rapid growth event that is not related to the skeletal core, and (2) this zone is not resolved well enough in the maps or the LA-ICP-MS data to be used as starting conditions for the modeling (i.e., exact shape and number of oscillations, P content of each oscillation). If this outer P-rich region is in fact related to the high-P skeleton, including it in the numerical model would improve the fit to the data but would not change the retrieved timescale. Measured LA-ICP-MS profiles were compared to a zone of similar dimensions and location on the modeled map: The location of the profile was estimated visually, and the Ti and Sc excesses along the profile were integrated over a 100-μm-wide zone reproducing the LA-ICP-MS acquisition method. Due to scatter in the LA-ICP-MS zoning profiles, the model’s fit to the data was estimated visually and not computed at each time step. The Fourier number, Fou, corresponding to the best fit was converted to a timescale using appropriate values for the diffusivity.

Fig. 8 2D lattice Boltzmann modeling of Ti and Sc zoning in the core of crystals 1957-A3 and B1. a-b Initial conditions, best fit profiles, and compositional maps for crystals A3 and B1, respectively. c Time elapsed between skeletal crystal growth and eruption as a function of temperature. The most likely range in temperatures and their associated timescales are shown with gray shaded areas. The temperature scale on the horizontal axis stretches from the coolest (905 °C, hydrous) to the warmest (1175 °C, anhydrous) temperatures measured for the three eruptions Full size image

The temperature-dependent diffusivity of Ti in olivine was studied by Cherniak and Liang (2014) under anhydrous conditions (at 1 atm with no addition of hydrous species). They observed no dependence on crystallographic orientation, Fo content, or oxygen fugacity and found very low diffusivities that are close to those of P (Spandler and O’Neill 2010) or Si in anhydrous olivine (Dohmen et al. 2002; Fei et al. 2012). Ti is in tetrahedral coordination and occupies the Si site in anhydrous olivine, but shifts to octahedral coordination and occupies an M1 site sharing an edge with a vacant Si site in hydrous olivine (Berry et al. 2007). Therefore, one would expect the diffusivity of Ti to be affected by that of Si in all cases. As the presence of only 10 ppm H 2 O in olivine will increase the diffusivity of Si by approximately three orders of magnitude (Costa and Chakraborty 2008), it seems likely that the diffusivity of Ti would increase by similar amounts in the presence of water. This would be the case for Llaima olivines, which host magmas contain ~2 wt% H 2 O. Crystals 1957-A3 and B1 display a diffusion length for Ti that is shorter than that of Fe–Mg at the rim and larger than that of P in the core (Fig. 6), in agreement with intermediate diffusivities that are greater than the experimentally derived anhydrous ones due to the stabilization of defects by OH groups (Berry et al. 2007; Costa and Chakraborty 2008). The diffusivity of Sc in olivine is not well constrained and was therefore taken to be the same as that of Ti.

Two-dimensional diffusion modeling cannot reproduce the measured profiles when initial enrichments dictated by the measured P enrichments are considered (i.e., Ti and Sc enrichments of 10.4 and 4.4 for crystal 1957-B1 and 4.8 and 3.2 for crystal 1957-A3). This could be due to the fact that Ti and Sc cannot be increased by such large amounts in the boundary layer melt due to structural considerations, or that diffusion in the third dimension has accounted for part of the decrease in concentration. 2D diffusion modeling can reproduce the Ti and Sc zoning profiles fairly well when lower enrichments of 4.5 and 2 for crystal 1957-B1 and 2.14 and 1.5 for crystal 1957-A3 are used (Fig. 8). Timescales for Sc are systematically larger, suggesting that its diffusivity is probably greater than that of Ti. Despite the fact that the timescales obtained here are very loosely constrained (scatter in the data—particularly Ti in crystal 1957-A3—and uncertain diffusivities and initial concentrations), we believe that the method has a lot of potential to determine older timescales from crystal interiors. An improvement of the current knowledge on diffusivities of trace elements such as Ti, Sc, or V in olivine, and the analysis of such elements with higher precision (e.g., by secondary ion mass spectrometry) would yield much more reliable timescales.

The skeletal core of crystal 1957-B1 likely formed within the recharge magma, hence its Fo content of 81 mol%. Temperature estimates for the core of this crystal are at least 1070 °C, considering hydrous conditions (1.5 wt% H 2 O measured in the melt inclusion). At this temperature, the skeletal core would have formed ~11 years prior to the 1957 eruption. It is difficult to estimate a single temperature for the diffusion of Ti and Sc out of the skeletal core of crystal 1957-A3, as this crystal has experienced a more prolonged and complex history of magma recharge, growth, and resorption. As for crystal 1957-B1, the skeletal core likely formed within the recharge magma, at temperatures around 1070 °C or higher. But the temperature then fluctuated and decreased through time. Temperature estimates for the core (Fo = 73.5 mol%) of crystal 1957-A3 are around 1035 °C, considering hydrous conditions (2 wt% H 2 O measured in the melt inclusion). At temperatures intermediate between 1035 and 1070 °C, this skeletal core would have formed ~11–19 years prior to the 1957 eruption.

Olivine zoning, magmatic processes, and timescales at Llaima

In summary, individual olivine phenocrysts appeared to record only one magma recharge event if only Fo, Mn, Ca, and Ni zoning are considered, but may in fact have experienced a complex history, including multiple recharge events given their P, Sc, and Ti zoning. This complex history involving multiple events of crystal growth, dissolution, and re-equilibration is similar to what is seen in the textures and An contents of plagioclase phenocrysts. It suggests that these minerals are coeval and have similar residence times in the shallow conduit/reservoir system at Llaima. In this view, the plagioclase and olivine records are reconciled. The initial dichotomy of simple versus complex magmatic processes experienced by both minerals no longer holds true.

Fo zoning at the olivine rims suggests that mixing with a recharge magma occurs on the order of a few months to a few years prior to eruption, for the small number of crystals we have examined. Many olivine crystals from some historic lava eruptions record large and abrupt reversely zoned profiles, and some of these extend to the grain edge without a late normally zoned outer rim. Rim compositions converge on similar values in some lava thin sections, but in others a wide range of rim compositions and a great diversity of zoning patterns are documented (see Fig. 5 of Bouvet de Maisonneuve et al. 2013). These observations can only be explained if crystals that had resided in crystal mush came into contact with recharge magmas shortly before eruption.

Crystal residence times in the shallow reservoirs are on the order of a few years to decades, as estimated by Ti and Sc zoning in their cores. These timescales suggest that the frequency of recharge events is on the order of a decade or more. This is less than the volcano’s overall frequency of eruption ~5.6 years (Dzierma and Wehrmann 2010), but individual olivines in the subvolcanic conduit/reservoir system probably record only a small fraction of the total number of recharge events, which must vary in magnitude and are probably dispersed among a large number of dikes and sills. Each eruption is probably associated with one or a swarm of variable-volume magma recharge events, in agreement with the large number of reversely zoned olivine phenocrysts in the tephra and lavas (Bouvet de Maisonneuve et al. 2012, 2013). Eruption volumes are likely to be crudely proportional to the aggregate volume of recharge magma in such swarms, and large eruptions may reflect the arrival of an anomalously large batch of recharge magma that may serve as the proximal eruption-triggering event. Closer links to pre-eruptive recharge and the nature of eruption-triggering processes require integration of petrologic data and model results derived from erupted samples with pre-eruptive and syn-eruptive geophysical data.

The 2008–2009 eruption is the only Llaima event for which there is precursory geophysical monitoring data. The precursory seismic activity that characterized the mid-2007 period had diminished to a lower level by late 2007, and there was no seismic swarm immediately before the onset of the main eruption on January 1, 2008 (Hugo Moreno, personal communication 2010). The methods developed in this study should be tested and refined by characterization and modeling of crystals acquired by real-time sampling of a subsequent eruption of Llaima that will benefit from ongoing, post-2008 geophysical monitoring and comprehensive observations during the eruption.

Unlike phosphorus, which records multiple recharge events, Fo, Mn, Ni, and Ca seem to record no more than one recharge event. This is surprising as Fe and Mg diffusion in olivine is sufficiently slow to record processes operating on the timescales of decades (Costa and Chakraborty 2004; Morgan et al. 2004; Costa and Dungan 2005; Morgan et al. 2006). The absence of information on early magma recharge events could be due to the fact that the signal from earlier recharge events was erased by diffusion and/or interferences from preexisting zoning, and is no longer resolvable analytically (Fig. 9). This is possible if variations in concentration at the boundary are of small amplitude (i.e., limited compositional contrast between the residing and incoming magmas) and last for a limited amount of time. The durations of interactions between the resident and replenishing magmas are presumed to be much shorter (weeks–months) than typical olivine residence times in crystal mushes during quiescence (years). Such offsets in the timescales of resetting compositions versus diffusive relaxation allow for any minor amplitude “spikes and troughs” in zoning profiles to be rapidly diminished. The unexpectedly long timescales obtained by kinetic modeling of Fo zoning could suggest that the modeled zoning profiles were affected by multiple overlapping recharge events of short duration, and are not exclusively the product of the latest event. If this is the case, normally zoned olivines that originated from the latest recharge magma are the most faithful recorders of the time elapsed between the latest recharge event and eruption.