DFENS: Diffusion Chronometry Using Finite Elements and Nested Sampling

In order to reconcile petrological and geophysical observations of magmatic processes in the temporal domain, the uncertainties in diffusion timescales need to be rigorously assessed. Here, we present a new diffusion chronometry method: Diffusion chronometry using Finite Elements and Nested Sampling (DFENS). This method combines a finite element numerical model with a nested sampling Bayesian inversion, meaning that uncertainties in the parameters contributing to diffusion timescale estimates can be obtained and that observations from multiple elements can be used to better constrain individual timescales. Uncertainties associated with diffusion timescales can be reduced by accounting for covariance in the uncertainty structure of diffusion parameters rather than assuming that they are independent of each other. We applied the DFENS method to the products of the Skuggafjöll eruption from the Bárðarbunga volcanic system in Iceland, which contains zoned macrocrysts of olivine and plagioclase that record a shared magmatic history. Olivine and plagioclase provide consistent pre‐eruptive mixing and mush disaggregation timescales of less than 1 year. The DFENS method goes some way toward improving our ability to rigorously address the uncertainties of diffusion timescales, but efforts still need to be made to understand other systematic sources of uncertainty such as crystal morphology, appropriate choice of diffusion coefficients, initial conditions, crystal growth, and the petrological context of diffusion timescales.

sion in spinel; Sr in plagioclase) can provide information of long-term magma storage times on the order of hundreds to thousands of years (Cooper & Kent, 2014;Mutch, Maclennan, Holland, et al., 2019;Zellmer et al., 1999;Zellmer et al., 2000), whilst faster diffusing species (e.g., Fe-Mg interdiffusion in olivine) can offer insight to processes operating days to weeks (Moore et al., 2014;Hartley et al., 2016;Lynn et al., 2017;Mutch, Maclennan, Shorttle, et al., 2019), or even minutes to hours (e.g., H + diffusion in olivine) before eruption (Barth et al., 2019;Newcombe et al., 2020). However, the value of diffusion timescales is diminished without proper petrological context and the rigorous consideration of underlying uncertainties. Indepth petrological characterization is required in order to determine whether the diffusion timescales can plausibly be linked to specific petrological processes, physical processes, and ultimately volcano monitoring data. Petrological observations are also required to test whether assumptions about initial conditions, boundary conditions and intensive parameters are appropriate.
Linking magmatic processes to geophysical observations through time requires a robust treatment of the uncertainties associated with diffusion timescales. The Arrhenius relationship between temperature and elemental diffusivity means that uncertainties in temperature play a dominant role in controlling error estimates. Many diffusion studies account for the uncertainties of the methods used to estimate temperature such as phase equilibria geothermobarometers (Ruprecht & Plank, 2013), however, the uncertainties in other intensive parameters that control diffusivity, as well as parameters in the diffusion coefficients themselves, are often not properly considered. Furthermore, the uncertainty structure associated with diffusion coefficients is correlated (Costa & Morgan, 2010). Here, we present a Bayesian inversion method, known as Diffusion chronometry using Finite Elements and Nested Sampling (DFENS) for modeling diffusion of multiple elements for timescale estimation. DFENS combines a finite element numerical diffusion model with a nested sampling Bayesian inversion scheme. This can simultaneously account for observations from multiple diffusing elements and produces more robust uncertainty estimates by taking account of the covariance in uncertainty structure of the underlying diffusion coefficients. The DFENS approach will help to improve our understanding of the variability of diffusion timescales in a single eruption that is a reflection of different growth, storage, and transport histories. Moreover, if we can better constrain the uncertainties on diffusion timescales of individual crystals, then it may be possible to disentangle temporal variations in natural crystal populations.
Few studies so far have considered diffusion in multiple mineral phases that record common magmatic histories, which can then be used to test the robustness of different mineral geospeedometers. In the plutonic record, Ca-in-olivine and Mg-in-plagioclase speedometers have shown consistent results when used to estimate the cooling rate of the lower oceanic crust (Faak & Gillis, 2016). However, in volcanic settings, complex crystal cargoes often make it difficult to compare different geospeedometers as different phases can record different magmatic histories (Chamberlain et al., 2014). The products of the Skuggafjöll eruption from the Bárðarbunga volcanic system, Iceland, contain macrocrysts of olivine and plagioclase that have been compositionally mapped in detail and appear to share a common history of long-term storage followed by rapid rim growth (Neave, Maclennan, Hartley, et al., 2014). Textural and microanalytical evidence indicates that these crystals provide a means of testing the consistency of olivine and plagioclase geospeedometers.

Multi-Element Diffusion Using the Finite Element Method
Diffusion chronometry relies on solving some variant of Fick's second law through time from a set of pre-defined initial conditions until the model matches the observed compositional data. In many silicate minerals, the diffusivity of the elements of interest are often spatially variable. For example, Fe-Mg interdiffusion, Ni and Mn diffusion in olivine depend on forsterite content (Chakraborty, 1997;Holzapfel et al., 2007;Petry et al., 2004;Spandler & O'Neill, 2010), whilst the diffusivities of trace elements in plagioclase (e.g., Mg, Sr, Ba, K) have been shown to depend on anorthite content (Cherniak & Watson, 1994;Van Orman et al., 2014). A spatially dependent version of Fick's second law (Equation 1) is therefore required to model diffusion for elements in silicate minerals that have a composition dependence (Costa & Morgan, 2010;Crank, 1979): where C is the concentration of the element of interest, D is the diffusion coefficient (diffusivity) and t is time. Diffusive coupling between different trace elements can also create additional complexity (Costa et al., 2003). In the case of trace element diffusion in plagioclase, forms of the diffusion equation that account for the chemical potential of the trace element component and coupling with the anorthite component need to be considered (Costa et al., 2003): where C is the concentration of the trace element of interest, X An is anorthite content (mole fraction), A is the dependence of the chemical potential of the trace element in plagioclase on the anorthite content, T is temperature (K) and R is the universal gas constant (kJ K −1 mol −1 ). The complex diffusive behavior in most silicate minerals, coupled with changing boundary conditions and diffusion coefficients imposed by continually changing intensive parameters in magmatic systems (pressure, P; temperature, T; oxygen fugacity, fO 2 etc.) makes it very difficult to solve diffusion timescale problems in igneous petrology using analytical solutions. This has led many studies to use numerical models to solve the diffusion equation using either finite differences (Costa et al., 2008;Druitt et al., 2012;Moore et al., 2014) or finite elements (Mutch, Maclennan, Holland, et al., 2019;Mutch, Maclennan, Shorttle, et al., 2019) that have been discretized in space and time.
The finite element method has emerged as a universal method for the solution of partial differential equations, like the diffusion equation. The power of the finite element method lies in its generality and flexibility allowing a wide range of partial differential equations to be solved within a common framework (Logg et al., 2012). A finite element is defined as a cell with a local function space (U) and rules that describe the functions that operate in this space (Brenner & Scott, 2008;Logg et al., 2012). Together these cells form a mesh which defines a functional domain (Ω). These meshes can take a range of simple polygonal shapes such as intervals, triangles, quadrilaterals, tetrahedra or hexahedra, which makes it a more useful way to generate complex morphologies such as crystal forms than regular finite difference methods ( Figure 1).
Here, we use the FEniCS software (Alnaes et al., 2015;Logg et al., 2012) to solve Equations 1 and 2. For this to happen, the unknown function (known as a trial function) needs to be discretized using the finite element method. This discretization involves multiplying the partial differential equation for the trial function by a test function (here represented as u) and integrating over the domain. Second-order derivatives are typically (but not always) integrated by parts. This new form is known as the "variational form" or "weak form" and holds for all u in some function space (U x ). The trial function (defined as C here for concentration) resides in a (possibly different) function space (U). These function spaces are defined by the mesh and the type of finite elements. A derivation of the variational form for a time-dependent diffusion problem is included in the Supplementary Material. The variational form for diffusion equations with a spatially dependent diffusion coefficient and time discretized according to a Crank-Nicholson scheme (Equation 1) is: where C k is the concentration at the previous time step k, C k+1 is the concentration at the next time step k + 1, C mid = (C k + C k+1 )/2, D(C mid ) is the compositionally dependent diffusion coefficient, Δt is the time step, u is the test function and Ω represents the spatial domain. The variational form used in this study for the plagioclase diffusion equation (Equation 2) is: where C mid , C k , C k+1 , Ω, Δt, u, R, X An , D, and A are defined above. For solving time-dependent partial differential equations the time derivative needs to be discretized by a finite difference approximation, which yields a recursive set of stationary problems that can then be written in variational form. We opted to use a Crank-Nicholson scheme because it is both stable and accurate. The trial function and the test function use the same functional space defined by the mesh and the type of finite element. A significant advantage of FEniCS is that it automatically does all of the discretization once the weak form has been characterized. This means models can be rapidly developed and are adaptable to complex problems. Once the partial differential equation has been discretized and finite element functional spaces have been assigned, the FEniCS software uses direct or iterative LU solvers to solve the resulting algebraic systems. For non-linear equations like Fe-Mg interchange in olivine, a Newton solver can be used. In all cases in this study, linear Lagrange (Continuous Galerkin) finite elements were used to represent concentrations.

Accounting for the Covariance in Uncertainty Structure in Diffusion Coefficients
Diffusion coefficient parameters are typically extracted using regressions through experimental data in ln D versus 1/T space via the Arrhenius relationship: MUTCH ET AL.
10.1029/2020GC009303 4 of 28  (Mutch, Maclennan, Holland, et al., 2019). (d) is a 2D finite element mesh of the crystal shown in (c) The mesh shown in (d) has been refined at its edges (i.e., has a smaller mesh size) so that a more detailed solution can be captured in areas of interest, such as where diffusion is most likely to be operating. This means a balance can be made between spatial resolution and computational time. BSE, back-scatter electron. where D 0 is the pre-exponential factor and E a is the activation energy. The slope and intercept of a linear regression are related to each other, which is critical when considering the uncertainties relating to the parameters that determine diffusion coefficients. This is particularly true for D 0 and E a , where higher values of D 0 would need to be associated with higher values of E a (Figure 2). Taking account of this form of uncertainty in diffusion modeling requires an understanding of the covariance of all the parameters that go into the diffusion coefficients. This feature has somewhat been neglected by most diffusion modeling studies. The main focus of this work is the creation of new multiple linear regressions through the experimental data so that the uncertainty structure can be properly assessed with covariance matrices. These regressions and covariance matrices are presented below and in the Supplementary Material, along with new modeling methods that can account for the trade-offs between different parameters.
New multiple linear regressions through a compiled database of olivine diffusion experiments (Chakraborty, 1997;Holzapfel et al., 2007;Petry et al., 2004;Spandler & O'Neill, 2010)  via the TaMED mechanism in olivine. The data were filtered for pressures at 1 atm, log 10 fO 2 at −7 Pa, and X Fo between 0.88 and 0.92. The inset is a density plot showing the covariance between these two parameters. A steeper gradient (−E a ) will be associated with a higher intercept (log 10 D 0 ), which is an important factor to consider for error propagation. For this example, the effects of olivine composition, pressure and oxygen fugacity have not been considered, but they are considered in the multiple linear regression presented in Equation 6.
where Ol, is the diffusion coefficient of species i in olivine parallel to the [001] direction, whilst a i , b i , c i , q i , h i and j i are the best fit parameters from the regression. X Fo is the forsterite content of the olivine (mole fraction). Pressure (P) is expressed in Pa, T in K and ln fO 2 in its native form (i.e., fO 2 is in bars). Versions of these equations with fewer parameters (i.e., no h i parameter is included) are also available in the Supplementary Material. It is important to note that the parameters shown in Equation 6 may be different to those that have traditionally been extracted from diffusion experimental studies (e.g., D 0 and E a ) as this study's regressions aim to fit all parameters simultaneously, whilst experimental studies often treat each parameter individually. Our regressions explicitly account for compositional effects (e.g., fO 2 and X Fo ) that are often wrapped up in the pre-exponential factor D 0 . Diffusive anisotropy is taken to be six times faster along the [001] axis than the [010] and [100] axes for Fe-Mg and Mn (Chakraborty, 2010), and 10.7 times faster for Ni (Spandler & O'Neill, 2010). In this study, we do not account for any uncertainties in diffusive anisotropy.
The covariance matrices associated with the fitting parameters from these new regressions are shown in the Supplementary Material. They were created so that the uncertainty structure associated with the experimental fits can be rigorously explored. As Mn is believed to diffuse via a similar mechanism to Fe-Mg interdiffusion (Chakraborty, 2010), Fe-Mg diffusion experimental data were used to supplement Mn data in order to determine Mn's diffusive dependence on forsterite content. The regressions recover all of the experimental data within 0.5 log 10 units and are consistent with previously reported equations (Chakraborty, 2010;Costa & Morgan, 2010;. The regressions and covariance matrices for Fe-Mg interdiffusion only use data from anhydrous experiments, and do not account for the effect of water on diffusivity (Hier-Majumder et al., 2005). The experimental data for Fe-bearing olivines show that activity of silica ( SiO 2 a ) only plays a minor role in Fe-Mg interdiffusion , and experiments for Ni and Mn have not been explicitly buffered for SiO 2 a so its effect is currently unknown. Separate regressions and covariance matrices for diffusion of Ni and Mn along [001] in pure forsterite from experimental data sets that were explicitly buffered for the activity of silica (Jollands et al., 2016;Zhukova et al., 2014) are included in the Supplementary material. The multivariate linear regressions performed for trace element (e.g., Mg, Sr, Ba, K) diffusion in plagioclase are presented using the form: where the regression parameters a i , b i , c i , and q i for the diffusion coefficient of species i in plagioclase ( Pl i D ) are not the same as those presented in Equation 6. These constants are presented in the Supplementary Material. The diffusion coefficients of Sr, Ba and K in plagioclase have dependences on anorthite content (Cherniak, 2002;Cherniak & Watson, 1994;Giletti & Casserly, 1994;Giletti & Shanahan, 1997), whilst the role of SiO 2 a has yet to be properly explored. For these elements c i would be equal to zero.
The compositional dependence of Mg diffusion in plagioclase has been explored in two experimental studies. The data set of Van Orman et al. (2014) considers the effect of anorthite content on diffusivity, but not the role of but does not include a term for anorthite content. We have also created an additional multiple linear regression through both data sets in an attempt to account for the effects of both of these compositional variables. Van Orman et al. (2014) report that all of their experiments were conducted under silica-saturated conditions, meaning we could assign them an SiO 2 a of 1 and that they can be potentially combined with the data of Faak et al. (2013).
Diffusive anisotropy has been shown to only play a minor role for most trace elements. For Mg it is thought to be approximately a factor of 2 (Van Orman et al., 2014), whilst no anisotropy has been reported for Sr (Cherniak & Watson, 1994). Our regressions include all data regardless of crystallographic direction and do not account for any of the effects of anisotropy between the [010] and [001] directions.

Parameter Estimation and Exploring Uncertainty Structure Using Bayesian Inference
We use Bayesian inference to directly estimate the parameters that contribute to our understanding of magmatic timescales based on multi-element diffusion chronometry. This method also provides a powerful way to explore the underlying uncertainty structure and for comparing the statistical likelihood of different physical models.
Bayesian inference is a method of statistical inference in which Bayes' theorem is used to update the probability for a hypothesis (or model) as more information, or evidence, becomes available. It involves calculating a posterior probability (the probability of a hypothesis given the evidence) from a prior probability (the probability of the hypothesis before the evidence is observed) and a likelihood function based on a statistical model of the observed data. Bayes' theorem for model selection states (Feroz et al., 2009): where H w is one hypothesis, or model, out of w competing hypotheses whose probability may be affected by the data () and the set of parameters (θ). For diffusion chronometry purposes, the hypothesis or model describes the proposed behavior of the system. It could relate to the diffusion mechanism of the element of interest or the magmatic phenomena generating the diffusion profiles which could manifest in initial conditions, boundary conditions or intensive parameters. The data () is what has been measured or observed, which would correspond to compositional profiles measured across minerals or melts. The parameters (θ) that describe the model such as time, intensive parameters and diffusion coefficients are being inverted for. P(θ | H w ) ≡ π(θ) is the prior probability of the hypothesis (H w ) before the evidence is observed. It corresponds to the probability distributions of the model parameters before they have been assessed relative to is the posterior distribution, which is the probability distribution of the parameters given the data and the competing hypotheses or models.
 is called the likelihood; it indicates the compatibility of the evidence with the given hypothesis. In this study, we define the following log-likelihood function: where obs  is the measured value, calc  is the value predicted by the forward model, and   is the standard  is the Bayesian evidence, which is the factor required to normalize the posterior over θ (Feroz et al., 2009): where N is the dimensionality of the parameter space. The Bayesian evidence inherently implements Occam's razor so that a simpler theory with a more compact parameter space will have a larger evidence than a more complicated one, unless the latter is better at explaining the data (Feroz et al., 2009).
The prior distributions can be described using different functions; the main ones used in this study are log uniform priors, Gaussian priors and multivariate Gaussian priors. A uniform prior is a constant probability function, which means that all possible values are equally likely a priori. A log uniform prior is a uniform prior that is applied across a logarithmic domain. In the models used in this study, time was assigned a log uniform prior due to the exponential relationship between temperature and diffusivity. A Gaussian prior uses a Gaussian probability distribution as defined by the mean and standard deviation. Intensive parameters that have been independently estimated, such as temperature (T), pressure (P), ferric iron content of the melt (Fe 3+ /Fe total ) and the activity of silica ( SiO 2 a ), were assigned Gaussian priors using the independent estimate as the mean and the inherent uncertainty of the method as the standard deviation. It should also be noted that thermobarometric methods may also introduce correlation between intensive parameters. A multivariate Gaussian prior involves the generalization of one dimensional Gaussian priors up to higher dimensions. This can account for any covariance in parameters (described by covariance matrices), which is the case for the parameters that contribute toward the diffusion coefficients. A series of univariate Gaussians can be converted into a multivariate Gaussian using: where m is the multivariate Gaussian, λ is a diagonal matrix of the eigenvalues of the covariance matrix, ϕ is the matrix of eigenvectors from the covariance matrix, ω is a one dimensional standard Gaussian distribution and μ is a vector of the mean values of the Gaussian distributions. Figure 3 shows how the prior distributions for a synthetic diffusion profile are related to the posterior distributions. Using a Bayesian approach to diffusion modeling allows for observations from multiple elements in single or multiple phases to be considered simultaneously. Considering the covariance of all of the parameters in the diffusion coefficients offers a more robust way of accounting for uncertainties. This is critical when trying to reconcile geophysical and petrological observations in the temporal domain.

Nested Sampling and the MultiNest Algorithm
Our approach aims to combine parameter estimation with parameter error propagation by assessing the posterior distributions in the region of maximum likelihood, that is, where the model best fits the data. To do this, we apply a Monte Carlo Bayesian inversion where all of the parameters are being estimated. Here, we use the MultiNest algorithm (Feroz et al., 2009(Feroz et al., , 2013Buchner et al., 2014) which employs ellipsoidal nested sampling, a type of Monte Carlo algorithm in which a fixed size of parameter vectors or "livepoints" are sorted by their likelihood (Skilling, 2004). The number of livepoints are typically set to 400 in order to balance efficiency and accuracy. A large number of forward models are run using the finite element diffusion models, and their likelihood is assessed by the log-likelihood function (Equation 9). In each forward model, the parameters that are contained in the livepoints are randomly drawn from the prior distribution and are clustered into multi-dimensional ellipses. This form of clustering allows MultiNest to follow local maxima with ease meaning the parameter space can be efficiently explored, which reduces the number of forward model runs required (Buchner et al., 2014;Feroz et al., 2009Feroz et al., , 2013. The algorithm keeps drawing new points until one is found with a higher likelihood than the least likely point which is then removed (Buchner et al., 2014), allowing the algorithm to scan from the least probable to most probable zones. The algorithm terminates once convergence of the marginal likelihood is attained (i.e., Bayesian evidence), and the maximum likelihood is adequately characterized.
We implement MultiNest version 3.1 using the pyMultiNest wrapper (Buchner et al., 2014), which allows for efficient integration with the Python interface of FEniCS. The model was also programmed with a Message Passing Interface (MPI), for parallel computing on multiple nodes. The DFENS model currently requires high performance computing in order to complete models in a reasonable time. Supercomputer clusters would be required for more complex problems, such as using high resolution 3D meshes, to ensure convergence to a solution. As an example, a Lenovo Thinkstation with an Intel XEON microprocessor could complete 10,000 1D olivine simulations in under 20 min when using 30 cores.
Once generated, the posterior distributions provide important information on the model parameters and the underlying uncertainties. In well constrained problems, most of the posterior distributions correspond well with the prior distributions ( Figure 3). This indicates that the posteriors are being controlled by the priors, which is useful for error propagation. If the posterior distributions lie substantially inside of the prior distributions, then the posterior distributions are being controlled by the data. This is most notable for the time parameter, which is unsurprising given that diffusion is a time-limited process. If there is significant deviation away from the prior distributions, then this may provide useful information about how the physical or diffusion model needs to be changed.
In most cases, the median values of the parameter posterior distributions, notably time and temperature, were used for further analyses. The median parameters, however, may not necessarily be the same as the MUTCH ET AL.
10.1029/2020GC009303 9 of 28 combination of parameters that produces the maximum likelihood solution (maximizes ). The mean of the posterior distributions was not used because it may be influenced by outliers. Figure 4 shows the covariance between the a Fe-Mg term and q Fe-Mg term from Equation 6 for Fe-Mg interdiffusion in olivine. This highlights the importance of including covariance into error propagation as it can reduce the size of the parameter space that is being explored. Accounting for covariance in diffusion parameters can significantly improve the uncertainty estimates, which will depend on the mineral phase, diffusing elements and timescales of interest. For Fe-Mg interdiffusion in olivine and for magmatic timescales on the order of 1 year, the 1σ uncertainties can be reduced by a factor of 1.5-3 ( Figure 4). The improvements in the robustness of uncertainty estimates mean that diffusion timescales can be compared to other observations (e.g., geophysical observations) in the time domain with more confidence. (a) shows the posterior timescales distributions (kernel density estimates) for different olivine Bayesian inversion models using the DFENS method that were used to fit synthetic olivine profiles. The profiles were made parallel to [100] using a time of 300 days, temperature of 1190°C, Fe 3+ /Fe total of 0.15 and pressure of 0.35 GPa, with additional noise added based on typical uncertainties from EPMA conditions used in this study (X Fo ∼ 0.01, Mn ∼ 36 ppm, Ni ∼ 36 ppm). The gray line marks 300 days, which was used to produce the data. The red curve is a Fe-Mg TaMED diffusion model that assumes that the parameters that control the diffusion coefficient are independent. The blue curve is a Fe-Mg TaMED diffusion model that includes diffusion parameter covariance as defined by the covariance matrix shown in the Supplementary Material. The purple curve is a multi-element diffusion model (Fe-Mg, Ni, Mn) that also includes covariance structure. (b), (c), and (d) are multivariate kernel density estimations showing the trade-off between posterior distributions in a Fe-Mg (the intercept) and q Fe-Mg (1/T term) for Fe-Mg interdiffusion. These plots have been color-coded using the same scheme as in (a) It is clear that models that include a covariance structure between the diffusion parameters are associated with much smaller uncertainties. EPMA, Electron probe micro-analyzer. Sigmundsson et al., 2015). The Bárðarbunga-Veiðivötn volcanic system comprises an extensive set of fissure swarms that have propagated up to 115 km to the southwest and 55 km to the north-northeast of Bárðarbunga central volcano ( Figure 5). It is the second largest volcanic system in the Eastern Volcanic Zone (EVZ), and elevated magmatic fluxes have been associated with the putative center of the Iceland mantle plume (Gudmundsson & Högnadóttir, 2007;Jenkins et al., 2018). Within historical times alone, eruptions in the EVZ have accounted for approximately 82% (∼71 km 3 ) of the estimated eruptive volume on Iceland (Thordarson & Larsen, 2007). During this period of time the Bárðarbunga-Veiðivötn volcanic system erupted at least 24 times making it the second most active system in historical time and therefore an important target for hazard management (Caracciolo et al., 2020;Larsen, 2002). The most recent Bárðarbunga-Holuhraun eruption in 2014-2015 serves as an additional reminder of the active nature of this volcanic system and the regional hazards that it can pose (Ilyinskaya et al., 2017;Sigmundsson et al., 2015;Ágústsdóttir et al., 2016).

Application of DFENS to a Petrologically Well
Deep seismicty (12−25 km depth) was detected beneath Bárðarbunga volcano up to 4 years before the Holuhraun eruption (Hudson et al., 2017). In the lead up to the eruption itself there was 13 days of seismicity that progressively propagated northeast from Bárðarbunga along the Dyngjuháls fissure swarm, which has been interpreted to represent the lateral propagation of magma (Sigmundsson et al., 2015;Ágústsdóttir et al., 2016). The eruption was accompanied by gradual caldera collapse, which supported the notion of lateral magma migration from the central volcano (Gudmundsson et al., 2016). The excellent coverage of geophysical monitoring methods of the Holuhraun eruption has provided a valuable insight into the timescales and mechanisms of dyke propagation and lateral magma flow during an Icelandic rifting event (Woods et al., 2018;Ágústsdóttir et al., 2016). These geophysical observations are now starting to be reconciled with geochemical observations in order to place real-time observations into a petrological framework Halldórsson et al., 2018;Hartley et al., 2018). However, to develop effective forecasting strategies for volcanic eruptions and their associated hazards, studies into multiple eruptions from the same volcano or volcanic system are required. In this instance, looking for pre-eruptive signals prior to dyke propagation in the petrological record of older eruptions may help to focus current geophysical monitoring methods of Icelandic volcanoes.
The Bárðarbunga-Veiðivötn system is also believed to have been highly productive during the Holocene and Pleistocene with large fissure eruptions repeatedly taking place on the south-western Veiðivötn fissure swarm (Larsen, 1984). The Skuggafjöll eruption is one such example of Pleistocene activity in the Bárðarbunga-Veiðivötn system. Skuggafjöll is an 820 m high mountain that is part of a NE-SW striking hyaloclastite ridge situated between Vatnajökull and Mýrdalsjökull (Neave, Maclennan, Hartley, et al., 2014). It is composed of plagioclase ultraphyric basalts that transition from pillow lavas at the base to hyaloclastites halfway up the mountain. These characteristics indicate that Skuggafjöll was a subglacial eruption, and places a minimum eruption age of approximately 10 ka (Jakobsson & Gudmundsson, 2008;Neave, Maclennan, Hartley, et al., 2014). A minimum erupted volume of 0.2 km 3 was estimated for Skuggafjöll by Neave, Maclennan, Hartley, et al. (2014) assuming a cone shaped edifice with a basal radius of 1 km and height of 0.2 km; although this did not take into account any subsequent erosion or burial by later eruptions. In spite of the poor constraints on eruption age and erupted volume, the well constrained petrological history preserved in its crystal cargo can be used to gain important constraints on the timescales of pre-eruptive processes in the Bárðarbunga system and to test the performance of different mineral geospeedometers.
MUTCH ET AL.
10.1029/2020GC009303 11 of 28  (EVZ) showing the location of the Skuggafjöll eruption (black diamond) within the Bárðarbunga-Veiðivötn volcanic system. The most recent eruption in the Bárðarbunga system, the 2014-2015 Holuhraun eruption, is also shown in purple for reference. The dyke propagation pathways for each eruption are shown as red arrows. For Holuhraun the dyke propagation pathway was constrained using pre-eruptive seismicity (Sigmundsson et al., 2015;Ágústsdóttir et al., 2016), whilst for Skuggafjöll a simple linear dyke pathway was assumed. The location of major central volcanoes is marked with their associated calderas (dashed lines). Major fissure swarms in the EVZ are shown in red (Thordarson & Larsen, 2007). Inset shows the location of the mapped region and Skuggafjöll with respect to the rest of Iceland. EVZ, Eastern Volcanic Zone.

Petrology and Sample Description
All samples described by Neave, Maclennan, Hartley, et al. (2014) of the Skuggafjöll eruption are olivine (1%-3%), clinopyroxene (2%-9%), and plagioclase phyric (3%-36%) with macrocrysts of these phases occurring as single isolated crystals and within monomineralic and polymineralic glomerocrysts. Plagioclase and olivine are often intergrown in glomerocrysts with interstitial melt pockets, which is suggestive of sequestration in a crystal mush as opposed to being joined by synnuesis just before eruption. The habit of many of the coarser plagioclase macrocrysts is too equant to be the result of rapid crystallization and is likely to represent a deep mush origin (Holness, 2014).
Whole rock geochemical variation indicates significant crystal addition, particularly of plagioclase (Neave, Maclennan, Hartley, et al., 2014). Olivine macrocrysts range in size from 150 µm up to 4 mm, and are typically equant and subhedral. Clinopyroxene macrocrysts are 150 µm to 2.2 mm in size with equant and prismatic habits. The plagioclase macrocrysts show the largest range in observed crystal size and texture. They range in size from 150 µm up to 12 mm with large, low aspect ratio (>600 µm size and length/width aspect ratios of 1.5) and small, high aspect ratio (<600 µm and aspect ratios >2) crystal populations present (Neave, Maclennan, Hartley, et al., 2014). Large plagioclase macrocryst cores show a range of melt inclusion textures from the absence of melt inclusions up to well-developed sieve textures. The presence of these defined crystal populations has been confirmed by crystal size distributions for each of the macrocryst phases, all of which show pronounced changes in gradient (Neave, Maclennan, Hartley, et al., 2014). The two crystal populations are also compositionally distinct; particularly for the cases of olivine and plagioclase. The coarser plagioclase and olivine macrocrysts have a more primitive character with core compositions of An 80−90 and Fo 85−87 , respectively. These crystal cores are surrounded by sharp, more evolved rims, An 70−79 and Fo 78−82 , that coincide with the compositions of the smaller macrocrysts and are in equilibrium with the matrix glass (Neave, Maclennan, Hartley, et al., 2014).
Melt inclusions from the primitive olivine and plagioclase macrocrysts show significant variation in their trace element compositions which is suggestive of crystallization from a suite of unmixed primary mantle melts (Maclennan, 2008;Neave et al., 2013;Winpenny & Maclennan, 2011). However, the major element composition of these different melt inclusion suites combined with the fact that their average trace element compositions are near identical within uncertainty provides strong evidence to suggest that the olivine and plagioclase cores co-crystallized from the same range of primitive melts (Neave, Maclennan, Hartley, et al., 2014). The average incompatible trace element composition of the melt inclusions is also significantly more depleted than that of the matrix glass, which indicates that the crystal cores and the more evolved rims crystallized from distinct melt distributions (Neave, Maclennan, Hartley, et al., 2014). Clinopyroxene-liquid geobarometry based on equilibria between the matrix glass and the clinopyroxene macrocrysts suggest that most crystallization took place at mid-crustal pressures (0.35 ± 0.14 GPa or 11 ± 4 km depth) (Neave & Putirka, 2017).
All of the above observations have been interpreted by Neave, Maclennan, Hartley, et al. (2014) to be the result of two stages of crystallization. The primitive macrocrysts cores crystallized from depleted primitive melts and were sequestered into a mineralogically stratified crystal mush pile in the mid-crust. Portions of non-cotectic mush were disaggregated and entrained into trace element enriched magma from which the more evolved rims and crystal assemblage grew at the three-phase gabbro eutectic. Transport and eruption at the surface must have occurred soon after given that the crystal rims are still relatively sharp. Modeling the diffusive re-equilibration between macrocryst cores and rims can provide a pre-eruptive timescale of the second stage of crystal growth and transport. The relatively simple petrological history that has been constrained by the in-depth work of Neave, Maclennan, Hartley, et al. (2014) makes Skuggafjöll an ideal eruption to develop, test and refine multi-element and multi-mineral diffusion modeling techniques.

Analytical Methods
Individual olivine and plagioclase crystals were picked from crushed glassy pillow basalt rims collected from the lower sections of the Skuggafjöll eruptive stratigraphy (GR: 63⋅968°N, 18⋅695°W). These were then mounted in epoxy 1-inch rounds and polished using silicon carbide papers and Metprep diamond suspension down to 0.25 µm grade.

Back-Scatter Electron Imaging
The texture and zoning patterns of approximately 40 olivine crystals and 50 plagioclase crystals were assessed by back-scatter electron (BSE) microscopy using a FEI Quanta 650FEG SEM at the University of Cambridge. BSE images were typically collected using an accelerating voltage of 10-20 kV and a working distance of 13 mm. To try to minimize charging effects from cracks and vesicles, 10 images were collected with a scanning rate of 1 µs and were integrated together with a drift correction. The brightness and contrast of collected images were adjusted using ImageJ image processing software in order to accentuate any zoning patterns. To minimize potential sectioning problems and diffusion from multiple dimensions (Costa & Morgan, 2010), crystal sections that followed the criteria of Shea et al. (2015) underwent quantitative analysis. Compositional profiles were positioned on euhedral crystal edges and in the center of crystal faces or as far away from other crystal edges as possible.

Electron Probe Micro-Analyzer
Compositional profiles of major and minor elements across selected olivine and plagioclase crystals were measured by electron probe microanalysis (wavelength dispersive X-ray spectroscopy, Electron probe micro-analyzer [EPMA]) using a Cameca SX100 with 5 wavelength dispersive spectrometers at the University of Cambridge. Calibration was carried using a mixture of natural and synthetic minerals and oxides. Instrument drift and measurement uncertainty was assessed by measuring secondary standards. For olivine analyses, an accelerating voltage of 20 kV was applied with a working current of 20 nA for major elements (Mg, Fe, Si) and 200 nA for minor and trace elements (Ni, Mn, Ca, Cr, Al). On peak count times of 20 s were used for major elements and 100-120 s for minor and trace elements, with half count times off peak. P was not measured routinely because the electron probe was operating without an LPET crystal (2 LIF arrangement). Plagioclase profiles were measured with an accelerating voltage of 15 kV and a working current of 10 nA for major (Ca, Al, Si, Na) and minor elements (Mg, Ti, K, Fe). On peak count times of 20 s were used for major elements and 90-110 s for minor and trace elements, with half count times off peak. For both sets of analyses, a spot size of 1 µm was selected, with profile point spacing varying from 5 µm (typically within 150 µm of the crystal edge) and 20 µm (distances exceeding 150 µm from the edge). For plagioclase, the beam was not defocussed to account for any alkali or silica drift given that Na and K concentrations were typically low in high anorthite plagioclase (Humphreys et al., 2006). Instead, Na and K were measured at the start of the analytical cycle for only 10 s.
Step scans (high resolution line scans) were collected by initially setting a line scan pre-sputter of 3.2 nA using 10 µm steps.
Step scan analyses were made with 2.5 × 10 −11 nA primary beam focused to approximately 2 µm, with step spacing set to 2 µm. There was no energy offset and 100 eV energy window was used. There were no losses due to field apertures as the spot size was much smaller than collection window. The scan position in the center of line was positioned with scanning ion imaging of Na and Si. Electron multiplier ions counting was used and all data were dead-time corrected (51 ns dead time). An entrance slit of 100 µm and exit slit of 400 µm were used. The nominal mass resolution was approximately M/ΔM 2400. A combination of feldspar (SHF-1 and Lake County plagioclase) and glass standards (NIST610, and V, W, X borosilicate glasses) were used to access analytical precision and convert raw counts to ppm values. Trace element silicon ratios measured by SIMS were then corrected relative to Si measured by EPMA.
Step scan data were then normalized to SIMS data in order to convert raw elemental ratios into concentrations. Prior to normalization, SIMS, step scan and EPMA profiles were projected onto a single profile that was orientated perpendicular to the edge of the crystal. Distances of analyses were corrected accordingly using cosΘ where Θ is the angle between the measured profile and the perpendicular projection.

Electron Backscatter Diffraction
Chemical diffusion of some major and minor elements in olivine has been shown to be strongly anisotropic. For example, Fe-Mg interdiffusion along the [001] direction is typically 6 times greater than along the [100] and [010] axes (Chakraborty, 2010;Costa & Morgan, 2010). The lattice orientations of the studied olivine crystals were thus characterized using electron back-scatter diffraction. Electron backscatter diffraction (EBSD) data with a resolution of 1-10 µm were collected at the University of Cambridge with a Bruker e Flash HR EBSD detector equipped on the Quanta 650FEG SEM, operating at 20 kV and beam spot size 5.5, and a stage tilt of 70°. The detector resolution was 320 × 240 pixels, while working distance and sample to detector distance were 17-30 and 12-18 mm, respectively. The data collection and indexing was performed with Bruker QUANTAX CrystaAlign software (QUANTAX, 2010), using a Hough transform resolution of 60-70. Data were analyzed using MTEX V4.0 (Bachmann et al., 2010), a freeware toolset for the commercial software package MATLAB (MATLAB, 2016).

Estimation of Intensive Parameters
The temperature of the carrier-liquid was estimated to be 1190 ± 30°C by Neave, Maclennan, Hartley, et al. (2014) using the clinopyroxene-liquid thermometer from equation 33 of Putirka (2008), which was applied to second generation clinopyroxene macrocrysts that were in equilibrium with the glass. A pressure of 0.35 ± 0.14 GPa was also estimated by Neave and Putirka (2017) using their recent clinopyroxene-liquid geobarometer. A Fe 3+ /Fe total (ferric iron content of the melt) value of 0.15 ± 0.02, representative of more enriched Icelandic basalts, was used (Shorttle et al., 2015); this value was then converted into an oxygen fugacity (fO 2 ) using an average glass composition of Neave, Maclennan, Hartley, et al. (2014) and Equation 7 of Kress and Carmichael (1991). These correspond to absolute ln fO 2 (bars) values of −18.76 ± 1 (QFM ± 0.3). The SiO 2 a (0.55 ± 0.04) of the Skuggafjöll magma was estimated using the same glass composition and the liquid's affinity for tridymite calculated in rhyolite-MELTSv1.02 (Ghiorso & Sack, 1995;Gualda et al., 2012).

Mg in Plagioclase Partitioning Behavior
Many of the empirical partitioning relationships (Bindeman et al., 1998;Nielsen et al., 2017) for Mg in plagioclase implicitly contain the dependence of the partition coefficient on temperature and melt composition in addition to anorthite content (Dohmen & Blundy, 2014). In order to try and isolate the dependence of the partition coefficient on anorthite content in the Skuggafjöll system, we adopt a similar approach as Moore et al. (2014) and focus on Skuggafjöll plagioclase macrocrysts with crystal faces defined by thin overgrowths. These rims are typically thinner than 20 µm (in some instances being only 5 µm thick) and are often associated with (010) faces that have slower growth rates than (001) and (100), respectively (Holness, 2014;Muncill & Lasaga, 1988). The parts of crystal cores adjacent to these rims likely equilibrated rapidly for Mg, meaning these faces provide an excellent opportunity to constrain the partitioning behavior of Mg in Skuggafjöll-like systems at a given temperature and melt composition. Rim and core compositional data measured within 20 µm of the crystal-melt interface were combined with experimental data (Bindeman & Davis, 2000;Bindeman et al., 1998) where K Mg is the partition coefficient of Mg in plagioclase and the numbers in brackets are the 1σ uncertainties on the fit parameters. The A Mg derived from this study has a negative slope, which is inconsistent with the thermodynamic analysis of plagioclase-melt partitioning data by Dohmen and Blundy (2014), which has a positive A Mg value. The nature of this discrepancy might depend on whether Mg preferentially partitions onto the M-site or tetrahedral site in calcic plagioclase (Dohmen et al., 2017;Longhi et al., 1976;Miller et al., 2006). Further work will be needed to understand the intricacies of Mg-in-plagioclase partitioning, however, for the purposes of this study, Equation 12 is suitable for application to Skuggafjöll.

Olivine Initial Conditions
Diffusion timescale estimates depend heavily on the assumed contributions of growth and diffusion, which is often expressed in the way that initial conditions are calculated. Compositional cross-plots of Al versus X Fo , Ni and Mn in Skuggafjöll olivines ( Figure 6) show step-like patterns that indicate potential diffusive decoupling between Al and the other diffusing elements. Experimental work by Zhukova et al. (2017) has shown that Al may diffuse rapidly via octahedral site vacancies, which is comparable to Fe-Mg interdiffusion. However, most Al in Fe-bearing magmatic olivine is incorporated in the tetrahedral site, and thus a slow diffusion mechanism coupled to Si is dominant (Spandler & O'Neill, 2010). Furthermore, in most of the profiles we measured, Al variation had a much shorter length scale than that of forsterite, and in some cases had sharp step-like morphologies (see Supplementary material). This suggests that the fast diffusion mechanism only played a minor role, and that the Al profiles are a feature of crystal growth rather than diffusion. Figure 6 also shows a convex pattern between X Fo and Ni, which indicates that most profiles were dominated by diffusion . Mutch, Maclennan, Shorttle, et al. (2019) assumed that Al profiles can be used to track the compositional morphology of rapid crystal growth and can thus be used as a proxy for initial conditions for the other elements of interest. This approach also relies on the assumption that the concentration of each element can be linearly related to each other during growth, and it is important to consider that this approach may not be applicable if zoning in Al and other elements are controlled by different processes. Textural and compositional observations made by Neave, Maclennan, Hartley, et al. (2014) show that the olivine rims crystallized concurrently with plagioclase and clinopyroxene following entrainment of crystal cores into the carrier liquid. We are therefore confident that for this eruption Al and Fe-Mg-Ni-Mn profiles in olivine are responding in a systematic way to this process, meaning we can MUTCH ET AL.

10.1029/2020GC009303
15 of 28 Figure 6. Compilation of olivine profile data collected by EPMA expressed as compositional cross-plots between the main elements typically used in olivine geospeedometry (X Fo , Ni and Mn) and Al, an immobile trace element (Spandler & O'Neill, 2010) that we use as a proxy for growth. The upper row corresponds to cross-plots between Al and X Fo (a), Ni (b), and Mn (c), whilst the lower row (d), (e) has Ni versus X Fo and Mn versus X Fo cross-plots. All of the data have been color-coded as a function of distance from the crystal edge. Cross-plots between Al and the elements of interest show a non-linear step-like distribution between rim and core compositions (purple lines) indicating diffusive decoupling. The large variability in Al content for forsteritic core compositions (X Fo ∼ 0.86-0.87) may reflect intercrystalline or intracrystalline heterogeneity in Al that has not been diffusively re-equilibrated in the crystal mush pile (Thomson & Maclennan, 2012). The cross-plot between Mn and X Fo shows a strong linear trend suggesting there has been very little diffusive decoupling between these two elements and that their diffusivities are similar. A subtle break in slope can be observed in the Ni versus X Fo cross-plot, which is indicative of minor diffusive decoupling likely imposed by slight differences in elemental diffusivity. Typical 1σ analytical uncertainties are shown by the black point. EPMA, Electron probe micro-analyzer.
adopt the same approach as Mutch, Maclennan, Shorttle, et al. (2019). Core and rim compositions of Al and the elements of interest were selected. Rim compositions were at the edge of the crystal and the core composition were chosen based on where the profiles flattened out (accounting for analytical uncertainties). A rim zone was selected based on where Al starts to decrease rapidly towards the edge of the crystal (taking into account any variations in Al content in the core). A linear calibration curve was then made between the rim and core compositions for each element. Diffusion would cause any deviations from linearity. The linear calibration curve was then used to convert Al compositions in the rim zone into concentrations of the element of interest. Points outside the rim zone were assigned the core composition. Figures illustrating this concept are in the Supplementary Material. As Phosphorus was not measured in most profiles, it was difficult to assess whether the Al profiles were controlled by growth rate. However, the fact that Al concentrations did not increase in the rim suggests that there was no enrichment associated with the establishment of a diffusive boundary layer (de Maisonneuve et al., 2016). Furthermore, the consistency in olivine rim compositions across all crystals (Al ∼ 160-180 ppm) suggests that rim composition may have been controlled by the far field melt composition.

Plagioclase Initial Conditions
Plagioclase initial conditions were developed using the assumption of the instantaneous growth of a rim in equilibrium with the surrounding melt. X An versus RT ln K Mg plots color coded for distance from the crystal edge (Figure 7) show that Mg compositions measured in plagioclase cores are negatively correlated with X An and form arrays that are subparallel to the partitioning relationship established in this study (Equation 12). Crystal rims and cores that are close to the rim-core interface typically fall off these trends which suggests that diffusion has taken place. These patterns indicate that the plagioclase cores equilibrated at a different set of P-T-X conditions (P-T-X 1) than those responsible for rim formation (P-T-X 2), with points between the P-T-X arrays representing disequilibrium. If the positive A Mg value of Dohmen and Figure 7. Calculated Mg partition coefficients (RT ln K Mg ) versus anorthite content for profiles collected by SIMS (squares) and EPMA (circles). Partition coefficients were calculated using the average concentration of the element in the glass and the estimated temperature of the carrier liquid (1190°C) (Neave, Maclennan, Hartley, et al., 2014). Each point is color-coded for the distance from the edge of the crystal. The gray lines are predictive partitioning models established for plagioclase at different sets of P-T-X conditions, that is, fixed liquid P, T, X (MgO content of the liquid) but variable X An content of plagioclase. The partitioning relationship is the one established in this study. The red arrow shows data that may have been influenced by diffusion. EPMA, Electron probe micro-analyzer; SIMS, Secondary Ion Mass Spectrometer. Blundy (2014) was applicable to Skuggafjöll, then the negative correlations in the core would need to be explained by plagioclase-dominated crystallization in which the effect of crystallizing mafic phases (e.g., olivine and clinopyroxene) on the melt Mg composition is negligible (Dohmen et al., 2017). Most MORB magmas, including Skuggafjöll, are expected to have crystallized along the plagioclase-olivine cotectic (in a 7:3 ratio of plagioclase to olivine), meaning that mafic phases still play a significant role in crystallization (Neave, Maclennan, Hartley, et al., 2014). Furthermore, Skuggafjöll plagioclase rims co-crystallized with olivine and clinopyroxene in eutectic proportions (Neave, Maclennan, Hartley, et al., 2014). This possibly rules out the role of crystallization in creating the observed negative dependence between anorthite content and Mg in the crystal cores. We interpret these signals to represent diffusive re-equilibration of plagioclase cores in a mush-like environment for a protracted period of time. This is supported by textural observations of mush storage and the homogenization of olivine compositions (Thomson & Maclennan, 2012). Mg initial conditions were produced by combining equilibrated core Mg compositions at P-T-X 1 conditions with a rim that was in equilibrium with the carrier liquid (i.e., there is a step in X An and the activity of Mg rather than continuous variation). The higher RT ln K Mg values calculated for core compositions suggest that they would be in equilibrium either at higher temperatures or with a more primitive melt (high MgO) than the final carrier liquid.

DFENS
Magmatic timescales were estimated for measured olivine and plagioclase compositional profiles using the DFENS method outlined above. A fixed Dirichlet boundary condition (C = C 0 on x = 0) was maintained at the crystal edge and a no-flux Neumann boundary condition (     0 on C x L n ) was maintained in the crystal interior. The standard number of mesh points for a 1D profile of length L was set to 300. The number of time steps in each realization was kept constant at 300; the size of the time step was not kept constant. The mesh was adapted and optimized according to the Courant-Friedrichs-Lewy (CFL) condition. Fe-Mg exchange was solved first at each time step using a Newton solver. Ni and Mn diffusion were then solved at each time step using the corresponding Fe-Mg (forsterite) solution. Diffusion of Mg in plagioclase was modeled using Equation 2. The models assumed that there was a semi-infinite melt reservoir.
A log uniform prior was used for time (10 −2 −10 4 days). Independent Gaussian priors, set with 1σ uncertainties, were used for intensive parameters including: temperature, pressure, ferric iron content of the melt, and the activity of SiO 2 . Multivariate Gaussian priors were used for coefficients in the diffusion equations that are controlled by their respective covariance matrices. In the case of plagioclase, a multivariate Gaussian prior was also used to define the A and B parameters of the Mg partitioning relationship (Equation 12) that contributes to the diffusive flux. This was constrained using the covariance matrix of the regression shown in Equation 12. The nested sampling Bayesian inversion was set with 400 livepoints, and the algorithm terminated once convergence of the marginal likelihood was attained.

Olivine Timescales
A total of 29 different olivine crystals were modeled using the DFENS method (e.g., Figure 9). The inversion typically converged to short magmatic timescales with the median of all modeled olivine crystals being 146 days and 95% of all retrieved timescales being shorter than 368 days (Figure 8). Each crystal generally required 10,000 to 300,000 realizations in order to reach convergence. The median values for all of the realizations for each individual modeled crystal ranges from 56 to 323 days. All of the olivine models converged around similar temperature, pressure and fO 2 conditions and are within the Gaussian priors used by the Bayesian inversion.

Plagioclase Timescales
Of the 22 plagioclase crystals modeled, 3 were not included in the final assessment due to uncertainties surrounding initial conditions and sectioning effects. In most cases the models provided good fits to the data (e.g., Figure 10). The resultant timescale distributions calculated using the DFENS methodology are dependent on the diffusion coefficients that were used. Plagioclase timescales calculated using the diffusion data of Faak et al. (2013) show excellent consistency with the olivine timescales. Figure 8 shows that the timescale distributions for these two phases are almost identical. The estimated median timescale is 140 days with 95% of timescales being less than 422 days. Timescales calculated using the Van Orman et al. (2014) data largely overlap with the olivine timescales. The median timescale from this distribution is 90 days, which is shorter than that for olivine, and 95% of timescales using Van Orman et al. (2014) are less than 219 days. For the regression that combined the data of Faak et al. (2013) and Van Orman et al. (2014), there is minor overlap on the upper bound of the olivine timescales. The median timescale of this distribution is 633 days whilst 95% of timescales are less than 2,118 days. The other intensive parameters, notably temperature, did vary more than those for olivine for each of the diffusion coefficients that were modeled. In some instances, they did converge outside of the original prior values. A Mg values ranged from −22 to −45; no models converged to positive values as suggested by Dohmen and Blundy (2014). Plagioclase crystals that converged to higher temperatures converged to lower A Mg values and vice versa. This could be due to the trade-offs between the trace element plagioclase partitioning relationships, which also controls the diffusive fluxes, and the other intensive parameters, most notably temperature.

Comparing Olivine and Plagioclase Timescale Estimates
Overall, there is good consistency between the timescale estimates obtained from olivine and plagioclase, particularly for plagioclase estimates using the separate diffusion coefficients of Faak et al. (2013) and Van Orman et al. (2014). Using the Mg in plagioclase diffusion coefficient that combines the data from both studies produces timescales that are typically four times longer than the olivine timescales. This discrepancy suggests that the data sets of Faak et al. (2013) and Van Orman et al. (2014) cannot be simply be combined. The two data sets likely form separate clusters that can be adequately described by individual linear regressions, however, a regression of the combined data sets has a significantly different slope. This could be due in their experiments. Even though free-silica was reported in the experimental charges, the SiO 2 a may not have been equal to 1 as we have assumed. Slight differences in diffusion mechanism could also account for discrepancies between experiments run at different anorthite contents. This complexity could relate to the sites in which most of the vacancy transport occurs (M-site vs. tetrahedral site) (Faak et al., 2013). Further study will be needed to reconcile differences between these two studies.
Given that the Faak et al. (2013) experiments were calibrated on plagioclases with anorthite contents (An 50−80 ) and bulk compositions ( SiO 2 a of 0.55-1) close to those observed in basaltic systems, we consider the diffusion data from Faak et al. (2013) as the best way for calculating Mg-in-plagioclase diffusion timescales in natural basaltic systems. We, therefore, base our interpretations on the timescales calculated using the Faak et al. (2013) data, which shows the greatest consistency with the olivine timescales. This suggests that rim growth took place less than a year prior to eruption. MUTCH ET AL.

Causes of Timescale Variability
The 1σ variation of both the olivine and plagioclase crystal populations is on the order of 200 days. Timescales for some individual crystals do not overlap within the uncertainty of the intensive parameters and diffusion coefficients calculated by the DFENS method. This variability could be the result of diffusion from multiple directions, sectioning effects, improper fitting, or uncertainties in partitioning models. These are discussed in more detail in the Supplementary material. Alternatively, the variation in timescales could be due to underlying magmatic processes.
Texturally, most olivine and plagioclase macrocrysts are very similar in that they have near homogeneous primitive cores surrounded by more evolved rims; this does make multiple magma storage regions unlikely, but does not preclude them. The plagioclase population does have subtle differences in trace element composition (e.g., Sr, Ba, La, and K) in their cores, but there is no relationship between core composition and pre-eruptive residence timescales. Some plagioclase macrocrysts that do have extra zones in their cores indicating that they experienced a more complex crystal history than that suggested by Neave, Maclennan, Hartley, et al. (2014). However, these crystals appear to have similar entrainment times to crystals with homogeneous cores.
Injection of new magma has often been invoked as a mechanism for initiating mixing and convection (Bergantz et al., 2015). Typical crystal residence times in the open convecting magma can be calculated following the method of Martin and Nokes (1989). This involves calculating a settling velocity for a spherical particle using Stokes' law: where v s is the settling velocity, Δρ is the density contrast between the crystal and melt, g is gravitational acceleration, α is crystal diameter, ρ is melt density and v k is the kinematic velocity of the melt. The settling velocity can then be combined with an exponential decay scheme to estimate the residence time (t r ): where τ is the thickness of the magma body. For a 10 m sill, a 2 mm diameter primitive plagioclase crystal (An 89 ) with a density of 2,641 kg m −3 would have a residence time of 160 days in a melt with a density of 2,704 kg m −3 and a kinematic velocity of 0.1 m s −1 . A 1 mm diameter primitive olivine crystal (Fo 86 ) of 3,285 kg m −3 density would have a residence time of 70 days. Crystal and melt densities are from Neave, Maclennan, Hartley, et al. (2014), which were calculated at 1190°C. The kinematic velocity was the upper limit for basaltic magmas from Martin and Nokes (1989). For a 100 m sill, the residence times for the same plagioclase and olivine crystals would be 1,500 and 700 days. It therefore seems that residence in a 10 m sill would be sufficient to account for the median diffusion timescales observed, though thicker magma bodies (∼100 m) would potentially be required to account for longer plagioclase residence times calculated via the combined diffusion equation. Additional complexity may arise from the fact that in some instances plagioclase and olivine cores are touching, meaning that there may be hindered settling or that some crystal clots are close to neutral buoyancy.
Incremental entrainment of crystal mush into the carrier liquid has been proposed as one mechanism for causing a range of observed timescales in basaltic fissure eruptions (Mutch, Maclennan, Shorttle, et al., 2019). This requires that the macrocrysts remain in contact with the magma for different periods of time. The duration of the Skuggafjöll eruption is unknown, however, given that many basaltic fissure MUTCH ET AL.
10.1029/2020GC009303 20 of 28 eruptions occur over months (Thordarson & Larsen, 2007), then this is the timescale over which diffusion in the open liquid could have plausibly taken place. Alternatively, the Skuggafjöll eruption itself may have taken place at the end of a much longer period of eruptive activity, although this is difficult to determine. Recent work by Cheng et al. (2020) that combines timescale estimates from diffusion chronometry with fluid dynamical simulations of magma intruding into crystal mush has shown a wide distribution of timescales can be associated with a single intrusive event. Crystals positioned in different parts of the remobilized mush may evolve along different P-T-X trajectories at different times, which may make it difficult to retrieve consistent timescales if these different conditions are not know a priori. Cheng et al. (2020) suggest that any delay between initial intrusion and when a diffusive response is recorded in the crystal cargo diminishes for longer magmatic residence times.
Neave, Maclennan, Hartley, et al. (2014) suggested that the non-cotectic character of the Skuggafjöll erupted products may have been the result of a mineralogically stratified mush. Plagioclase crystals concentrated at the top of the mush may have been preferentially entrained into the carrier liquid leaving behind olivines at the base. If this were the case, we would expect slightly different timescale distributions for the plagioclase and olivine assuming there was perfect segregation between the different phases. The models of Bergantz et al. (2015) and Cheng et al. (2020) suggest that material at the base of the mush pile would be exposed to the new liquid earlier on in the intrusive event and would thus have longer timescales. If we consider the plagioclase timescales from the Faak et al. (2013) diffusion coefficient, then the similarity between the olivine and plagioclase timescale distributions suggests that the plagioclase and olivine crystals may have been sampled from similar parts of the mush pile. This may suggest that there was not perfect segregation of olivine and plagioclase via mechanisms such as hindered settling or synnuesis. Sampling a larger part of the crystal population, minimizing uncertainties associated with sectioning and model fits, and reconciling Mg-in-plagioclase diffusion coefficients may help to further tease apart natural variation in pre-eruptive residence timescales and resolve potential discrepancies between the timescale estimates of olivine and plagioclase.

Placing Diffusion Timescales Into a Petrogenetic Context
The pre-eruptive timescales estimated in this study can be placed into the context of at least two phases of crystallization from geochemically distinct magma batches as proposed by Neave, Maclennan, Hartley, et al. (2014) (Figure 11). Primitive plagioclase and olivine macrocryst cores co-crystallized from primitive depleted melts at mid-crustal pressures (∼11 km depth). Trace element variability in olivine-hosted melt inclusions suggests that magma mixing was taking place concurrently with crystallization. The morphology of olivine-plagioclase contacts in glomerocrysts suggests that these crystals were then sequestered in a crystal mush rather than being joined by synnuesis (Neave, Maclennan, Hartley, et al., 2014). Diffusive equilibration of Mg in plagioclase cores and forsterite in olivine crystal cores suggests that this storage must have lasted at least a few hundred years (Cooper et al., 2016;Mutch, Maclennan, Holland, et al., 2019;Thomson & Maclennan, 2012). Following this period of protracted mush storage and re-equilibration, the mush was then disturbed and disaggregated by a more evolved melt that had originally differentiated at depth. This injection event would have accompanied the second phase of crystallization, and may have efficiently mixed with the host primitive magma if injection was rapid (Bergantz et al., 2015). The efficient mixing between the two liquids and the mush liquid for a long period of time could explain why no mush liquid component is observed when crystal addition is accounted for in the composition of whole rock samples (Neave, Maclennan, Hartley, et al., 2014). The entrainment of this mush into a now well mixed magma that is slightly colder would have promoted the observed rapid rim growth. Our petrological observations and diffusion timescales suggest that crystal residence in this newly mixed magma and transport to the surface took place less than 1 year before eruption. This helps to rule out a second petrological scenario proposed by Neave, Maclennan, Hartley, et al. (2014) which suggested there was shallow storage of evolved melts prior to eruption. This scenario would have required an extra phase of crystal growth and additional zones that are not observed. Furthermore, the volatile contents of olivine-hosted melt inclusions are likely the result of decrepitation upon ascent rather than representing shallow entrapment pressures (Maclennan, 2017).

Comparison With the 2014-2015 Holuhraun Eruption and Implications for Hazard Management
The final crystal entrainment and transport of the Skuggafjöll magma took place approximately 50-300 days before the eruption. Seismicity detected prior to the Holuhraun eruption indicate that magma transport time took place over approximately 13 days. This is corroborated by diffusive hydration timescales of olivine-hosted melt inclusions which provide a minimum estimate of magma residence time of 1-12 days ). An in-depth diffusion chronometry study has yet to be published on magmatic zoning of Holuhraun macrocrysts so crystal entrainment and residence in the final magma prior to the initial dyke propagation event are still unknown.
It is unclear whether dyke propagation and magma migration prior to the Skuggafjöll eruption would occur over similar timescales to that of Holuhraun. The distance between Bárðarbunga central volcano and the Skuggafjöll eruption site is approximately 60 km when assuming a linear propagation pathway. This distance is approximately 1.5 times the dyke propagation distance of Holuhraun, suggesting the timescales for Skuggafjöll are likely to be similar. Sigmundsson et al. (2015) have suggested that underlying topography and its influence on gravitational potential energy can play a large role in controlling the orientation of the dyke. This is particularly prominent close to the central volcano where topographic load is high, whilst regional tectonic stress fields play a more important role at distal portions of the propagating dyke tip. As Skuggafjöll was erupted during the last glacial period, when MUTCH ET AL.
10.1029/2020GC009303 23 of 28 Figure 11. Schematic cartoon showing our proposed model for the petrogenesis of the Skuggafjöll magma, which involves two stages of crystallization. Olivine is shown in green and plagiolcase in white. (a) shows the crystallization of the primitive macrocryst assemblage from geochemical variable melts (first stage of crystallization). (b) shows the sequestration of these primitive macrocrysts in a crystal mush. The second stage of crystallization is outlined in (c) and (d). Recharge of the primitive mush with a more evolved and enriched magma (c), causes plagioclase dissolution and mush disaggregation, followed by the second stage of crystallization prior to eruption (d). Diffusion chronometry using DFENS suggests this second phase of crystallization and mixing took place approximately one year before eruption. Figure adapted from Neave, Maclennan, Hartley, et al. (2014). DFENS, Diffusion chronometry using Finite Elements and Nested Sampling.
there was additional loading of the crust by glacial ice, modern day topography may be ill-suited for predicting the dyke pathway leading to the eruption site. Regardless, any changes in dyke propagation path are likely to be minor as most of the pathway was distal from the central volcano and would thus be controlled by tectonic stresses, which is close to the down rift linear approximation. Any modification in transport time is therefore likely to come from the dyke stalling in the crust, which cannot be determined. Any lateral or vertical magma transport to Skuggafjöll is unlikely to have taken more than a few weeks, meaning most of the timescale recorded by the crystal cargo probably relates to mush reorganization and magma movement at depth.
Deeper seismicity (12-25 km depth) to the east of Bárðarbunga was detected up to 4 years before the Holuhraun eruption (Hudson et al., 2017), which could be interpreted as magma mixing and supply of melt from greater depths. The timescales and depths of this activity and that estimated from the crystal record of Skuggafjöll make for a tempting comparison given that they are fairly similar (i.e., deep activity recorded months before eruption). It could be speculated that that these events refer to a common process (i.e., melt migration from greater depth followed by magma mixing and crystallization), however, the lack of geophysical observations prior to Skuggafjöll and lack of diffusion studies of Holuhraun mean that a model of magma emplacement and mixing months to years before eruption would require more multi-disciplinary observations in order for it to be applicable for forecasting basaltic fissure eruptions.
A further note of caution for comparison relates to differences in melt inclusion trace element compositions between the two eruptions. The composition of olivine-hosted Skuggafjöll melt inclusions  is typically more depleted than that of Holuhraun and other eruptions from the Bárðarbunga system . This is in spite of the fact that the whole rock compositions fall within the Bárðarbunga-Veiðivötn array. This may suggest that Skuggafjöll was sourced from a slightly different part of the system. Nevertheless, if consistent deep pre-eruptive magmatic behavior can be shown for other case studies from the Bárðarbunga system, detecting deeper seismicity may be the strongest indicator that an eruption may be imminent within the following few years which may aide planning and hazard management in the area over this time period.

Conclusions
Diffusion chronometry applied to magmatic crystals plays a significant role in characterizing the temporal evolution of volcanic plumbing systems and reconciling geophysical and petrological observations. However, robust uncertainty propagation associated with this form of quantitative petrology has yet to be fully realized. A new Bayesian inversion method that combines a finite element numerical model with a nested sampling approach (DFENS) has been developed in order to achieve more robust uncertainty estimates and to account for the observations from more than one element within a single phase or multiple phases. This method offers a promising way to account for multi-element diffusion timescales from different minerals to be adopted into a single framework. We applied the DFENS method to olivine and plagioclase macrocrysts with a shared magmatic history from the Skuggafjöll eruption to estimate the timescale between crystal entrainment and eruption. There is excellent agreement between both phases which return timescales on the order of hundreds of days; olivine had a median time across all crystals of 146 days and plagioclase had a median of 140 days as calculated using the diffusion coefficient parameterization of Faak et al. (2013). The parameterization of Faak et al. (2013) may give the best timescale estimates for plagioclase residence because the data were calibrated at conditions closer to natural basaltic systems.
The estimated timescale of months to years for mush disaggregation and entrainment prior to the Skuggafjöll eruption could be comparable to deep seismicity detected up to 4 years before the 2014-2015 Holuhraun eruption, which has been interpreted as the upward migration of deep melts (Hudson et al., 2017). This highlights how the combination of detailed petrological work on erupted products, diffusion timescales with robust uncertainty estimates, and geophysical measurements of deep seismicity have significant potential in forecasting basaltic fissure eruptions.