Modeling of ionic diffusion in natural crystals has been developed over the last three decades to calculate timescales of geological processes. As the number of studies and the size of data sets have expanded, improvements in the precision of the general technique are needed to resolve temporal patterns that would otherwise be masked by large uncertainties. This contribution examines fundamental aspects of timescale calculation uncertainty using Mg-Fe zonation in olivine crystals from a Piton de la Fournaise oceanite erupted in 2002 CE. First, we quantitatively consider the role of geometric uncertainty in data sets from the perspectives of sectioning angle, crystal shape, and crystal agglomeration. Second, we assess how crystal growth and changing boundary conditions during diffusion pose problems for simplistic, 1D, diffusion-only modeling.

An initial database of 104 timescales (7–45 days) was generated using typical, 1D, isothermal diffusion-only methods for profiles taken from 30 compositionally and texturally zoned crystals of olivine. This simplistic modeling yields poor model fits and imprecise timescales; prior to this work, we would have rejected >60% of these data. Universal-stage measurements of crystal boundary angles and three-dimensional (3D) X-ray microcomputed tomography observations of crystal shape address geometric uncertainties. U-stage measurements show that, contrary to expectations of random sectioning, most boundaries modeled initially were close to the ideal sectioning plane. Assessment of crystal morphology from 2D thin sections suggests olivine crystals are dominantly euhedral; however, 3D imaging reveals that they are significantly subhedral and often exist as agglomerates, an observation that underscores both the potential for diverse crystal interactions through time in the magma (Wieser et al. 2019) and out-of-plane effects capable of influencing calculations of diffusion profiles.

Refinements to timescale determination can be made using a dynamic 1D modeling code to resolve growth and changing boundary effects simultaneous with diffusion. We incorporated temperature-dependent crystal growth rates (both linear growth and quadratically increasing, with a peak growth rate ~1.9 ×10–11 ms–1) and temperature-dependent boundary conditions (controlled using a cooling rate of −0.5 ± 0.1 °C/h) to remodel 13 timescales, resulting in significantly improved fits of the diffusion model to the initial data, better agreement between different faces of the same crystal, and less scatter within the whole data set. The use of 3D imaging and the inclusion of changing boundary conditions and crystal growth for diffusion calculations will enable more robust conclusions to be drawn from similar data in the future. Accurately retrieving timescale information from these crystals expands the pool of data available and reduces sampling bias toward “well-behaved” crystals.


Understanding of magmatic plumbing systems has been advanced in recent years by investigation of individual crystal histories using quantitative diffusion chronology (Costa and Dungan 2005; Morgan et al. 2006; Kahl et al. 2011; Charlier et al. 2012; Shea et al. 2015a; Hartley et al. 2016; Pankhurst et al. 2018b; Ruth et al. 2018; Petrone et al. 2018; Mutch et al. 2019). The essential observations underpinning these studies are that chemical disequilibria within natural crystals are often partially relaxed toward full equilibrium: a homogeneous distribution of elements (Fig. 1a). Timescales are calculated by modeling a diffusion profile that fits the observed profile (Costa and Morgan 2010).

Knowledge of ionic diffusivity in mineral lattices is gained from laboratory experiments (e.g., Buening and Buseck 1973; Chakraborty et al. 1994; Lesher 2010). These diffusivity measurements contain inherent uncertainty, which can be propagated to modeled results from natural data. Although experiments are designed to control environmental conditions and geometric relationships (Costa et al. 2008; Brady and Cherniak 2010; Chakraborty 2010; Dohmen and Milke 2010), it cannot be assumed that nature is so well-behaved, i.e., disequilibrium conditions may be prevalent across the system and subhedral, aggregated crystals may be present. Variations in intensive parameters such as temperature, pressure (Blundy and Cashman 2008), composition (Dohmen and Chakraborty 2007a, 2007b; Dohmen et al. 2007), and the geometric relationships of crystals (single crystals vs. aggregates, euhedral vs. subhedral) are central to how igneous rocks are formed (Welsch et al. 2013; Shea et al. 2015a). Application of experimental results to natural systems, while assuming these parameters are invariant, will likely lead to inaccuracies.

This contribution aims to improve the utility of diffusion chronology in natural rocks by investigating the effects of three common sources of analytical and epistemic uncertainty in natural systems: (1) sectioning angle, (2) 3D crystal shape, and (3) crystal growth and changing boundary conditions. We use Fe-Mg inter-diffusion in natural olivine as our case study, as it has a well-constrained and understood diffusion coefficient (Dohmen and Chakraborty 2007a, 2007b; Dohmen et al. 2007; Costa et al. 2008; Kahl et al. 2013; Bouvet de Maisonneuve et al. 2016; Lynn et al. 2017) and because basalt petrogenesis is relatively simple. The types of improvements we discuss, however, are equally applicable in other minerals and more complex igneous systems.

We constructed a data set of core-rim composition profiles using traditional methods (see Fig. 1b). From these profiles, a diffusion chronometry data set was generated whereby a single temperature is input into a 1D diffusion model, and solutions were corrected for anisotropy (Costa and Chakraborty 2004; Kahl et al. 2011; Hartley et al. 2016; Morgado et al. 2017; Pankhurst et al. 2018b) (Fig. 1c). We then made observations using a universal stage, X-ray microcomputed tomography (XMT), and a modeling program that can accommodate dynamic boundary and temperature conditions. Uncertainties within the chronometry data set are then defined and corrected by applying more advanced methods of observation and the advanced model.

We assess the influence of sectioning angle, crystal shape, and crystal growth, and changing boundary conditions at both a crystal scale and at a population scale. This work is complementary to the theoretical study of geometrical uncertainties and their effects on the accuracy of timescales demonstrated in Shea et al. (2015a). While Shea et al. (2015b) discriminated between growth and diffusion within a classic skeletal olivine crystal, further validation requires diffusion-based studies that aim to retrieve timescales from a wider pool of natural crystals. The studied olivines are at a late stage of crystal formation, the timing of the chemical profiles occurring after post-textural ripening to their observed form.

Sampling and geological setting

Piton de la Fournaise (PdlF) is a basaltic shield volcano forming the south-eastern portion of La Réunion, one of the Mascarene Islands that lie ~800 km off the coast of Madagascar (Longpré et al. 2006; Lénat et al. 2012) (see Fig. 2). There have been multiple eruptions in nearly all decades since historical records began in 1640 (Global Volcanism Program). Recent eruptions have occurred in 1990, 1991, 1992, 1998–2010, and 2014–present. The recent activity is mainly confined to the central caldera—Enclos Fouqué—and from the northeast and southeast rift zones (Lénat et al. 2012).

The November 2002 eruption (Fig. 2) occurred on the eastern flank as a lava flow (Longpré et al. 2006). The eruptive period is divided by Longpré et al. (2006) into five eruptive episodes defined by periods of inflation/deflation at the summit and the flanks and increased or decreased seismicity related to dike injection within the eastern flank (Longpré et al. 2006). The lava flow is comprised of oceanite [a term first introduced by Lacroix 1923 from Peltier et al. (2009)] and formally defined by Boivin and Bachèlery (2009).

We sampled the 2002 lava flow from a road-side outcrop at Lat. −21.13136S, Long. 55.48424E at 223 feet elevation. This lava contains abundant 2–3 mm olivine crystals. The sample collected is denoted “RU0701” throughout the paper, and whole-rock geochemical information for this sample, including major-element composition, was reported in Peters et al. (2016).

Initial methods

Sample preparation and characterization

Three thin sections were cut to 48 mm ×26 mm ×35 μm, ground, and polished with colloidal silica to give sufficiently defect-free finish suitable for electron-backscatter diffraction (EBSD) analysis, as in Lloyd (1987). Within each thin section (RU0701_1-3), we refer to an individual olivine crystal using a letter and use a number to refer to a unique compositional profile extracted from that crystal for the modeling, for example RU0701_1_A1. A scan of each thin section is included in Online Material1 Data A. Mineralogical modal abundance was calculated using “JMicroVision” as applied to scanning electron microscope images in addition to visual assessment of olivine shape made in both 2D and 3D. A core (diameter = 1.5 cm, length = 10 cm) was taken from the bulk-rock sample using a drill corer at the University of Leeds, in preparation for X-ray microcomputed tomography scanning.

Scanning electron microscope and electron probe microanalysis

Backscattered electron (BSE) images of thin sections were collected for textural assessment and for composition calibration using a FEI Quanta 650 Field Emission Gun-Environmental Scanning Electron Microscope (SEM) at the University of Leeds. The imaging was conducted using a focused spot size and an accelerating voltage of 20 kV. An Oxford Instruments electron backscatter diffraction system using a hkl Systems Nordlys detector was used to gather crystallographic orientation data from the portions of olivine used for compositional profiling. These data were reduced using Oxford Instruments Channel 5 software to a single set of Euler angle rotations per grain, and an in-house MS Excel worksheet was used to transform these into a set of a-, b-, and c-axis orientations relative to the EBSD analytical plane (i.e., the thin section plane).

Compositional rim-to-core traverses from the olivine grains were measured using a JEOL JXA8230 electron probe microanalyzer (EPMA) equipped with five wavelength-dispersive spectrometers (WDS) at the University of Leeds. A focused beam was used to extract a profile perpendicular to the crystal edge at 5–10 μm intervals. The broadest diffusion margin of the crystals measured is ~80 μm.

Multiple analyses were conducted to build the full data set, using slightly different beam conditions (15–20 kV and 30–50 nA) and on-peak count times (30–60 s for Ni and Mn; 30–40 s for Si, Fe, Al, and Mg; and 20–30 s for Cr). Higher energy and currents were used together with shorter count times and resulted in similar precision, as calculated using data from secondary standards from each run. Analyses that recorded totals outside 98 to 101.1 wt% were rejected. Microbeam reference materials distributed by the Smithsonian Institute, Washington, D.C. (Jarosewich et al. 1980) were used as primary standards, whose details along with detection limits determined in this study are reported in Online Material1 Data A.

The San Carlos olivine standard NMNH 111312-44 was used to calculate external accuracy and reproducibility of olivine composition via analyses at the start, end, and within the analytical session. The Fo content of the San Carlos olivine was reproducible within-analytical-session precisions of 2σ = 0.07–0.14 mol% (n = 66 over the course of the study). Internal olivine standards used were Springwater (USNM 2566) and Geo2 (see Pankhurst et al. 2016). Analyses of groundmass minerals were conducted alongside suitable matrix-matched standards, as detailed in Online Material1 Data A.

Initial modeling

Thermometry and oxygen fugacity.

For mafic magmas, it is common to utilize olivine-melt thermometry (Putirka 2008) to determine crystallization temperatures. In the absence of preserved homogeneous glass, we retrieved bulk compositions of the microcrystalline groundmass (Geiger et al. 2016). A quantitative mapping, raster matrix method was used to collect the compositional data, using the Probe for EPMA software (Donovan et al. 2012) and further refined using Surfer Software.

Grids of the groundmass were scanned using WDS, after which a matrix correction was applied to the data on a pixel-by-pixel basis, mimicking beam defocusing (Morgan and London 1996). The raster matrix method has the advantage of analyzing large areas without the beam attenuation and spectrometer viewing angle problems associated with using wide-area defocusing (Reed 2005). The MgO content (in wt%) derived from the raster matrix was input into the Helz and Thornber (1987) liquid-only thermometer. This thermometer has been widely implemented to investigate the thermometry of PdlF lavas in previous studies, e.g., Clague and Denlinger (1994); Famin et al. (2009); Bureau et al. (1998); and Boivin and Bachèlery (2009). Uncertainties presented in this study incorporate both the internal uncertainty of ±10 °C inherent to the thermometer as well as the MgO variability obtained via the groundmass analysis. Quantitative groundmass maps and calculation of their uncertainty can be found in Online Material1 Data A. We set the oxygen fugacity at the NNO+1 buffer. A detailed explanation of this choice is included in Online Material1 Data B.

Timescale calculations.

Profiles of grayscale values from BSE images were extracted perpendicular to the crystal edges using Image J (Schneider et al. 2012), illustrated in Figure 1b. The grayscale values were calibrated to the major element compositions (Costa and Morgan 2010) obtained from analyses of the same olivine crystal using the EPMA. Best-fit diffusion profiles were then modeled for these profiles using the same numerical methods as previous work by Hartley et al. (2016), Pankhurst et al. (2018b), and Couperthwaite et al. (2020). A demonstration and explanation of the methods of the simple diffusion model used are supplied in Online Material1 Data C.

Initial results

Olivine composition and textures in two dimensions

Whole-rock (XRF) compositional data for RU0701 (Online Material1 Data A), plots in the range for an ultra-mafic/basaltic (SiO2 wt% < 45%) rock with a high MgO content (23.5 wt%), consistent with the petrographic identification as oceanite [8–28 wt% MgO, Lacroix (1936)]. The large olivine crystals that are the main focus of this study appear as mostly euhedral, occasionally embayed, crystals up to 2–3 mm in size, at ~25% modal abundance. There are ~300–400 olivines in each thin section. They are contained within a finer-grained groundmass (~50–80 μm) made up of clinopyroxene (40%), olivine (10%), plagioclase feldspar (50%), and spinel (5%) (see Fig. 3), and often contain plagioclase inclusions in their margins. Olivine ~40–200 μm in size is also present (Fig. 3), much smaller than the population of large crystals.

The large olivines have core compositions of Fo84 [Fo = 100·(Mg/Mg+Fe+Ni+Mn)], consistent with previous olivine measurements for RU0701 (Füri et al. 2011), and Fo rim compositions ranging from Fo75–65 (normal zonation) with what is recognizable as diffused zones typically ~80 μm thick (Fig. 3). The olivine rims are similar to the groundmass olivine composition (~Fo65).

External modeling parameters and timescale uncertainties

A temperature of 1126 °C (±10 °C) was calculated using the average MgO (5.55 MgO wt%) of four quantitative groundmass maps (values range between 4.9 and 6.5 ± 0.5 wt% MgO; Online Material1 Data A). This estimated temperature is in good agreement with that used by Boivin and Bachèlery (2009), who calculated temperatures usually between 1110–1150 °C (and up to 1172 °C, with an interquartile range of 33 °C) for various eruptions from 1977 to 1998. We set the fO2 at NNO+1 log units.

The initial condition was defined using the core composition Fo84, with a variation in the boundary condition inferred from the rim composition between Fo75–65 depending on the profile being modeled. The average uncertainty of each calculated timescale is 0.38 log units (1σ), calculated considering the uncertainty on temperature, activation energy, and the pre-exponential factor that relates to crystal structure, jump frequency, and distance (Dohmen and Chakraborty 2007a, 2007b; Dohmen et al. 2007).

Initial Fe-Mg diffusion timescales

Diffusion timescales (n = 104) were calculated for Fe-Mg inter-diffusion in the large, zoned olivine crystals (n = 30). One profile was extracted from each available face of a single crystal and/or was taken along each diffusion direction (relative to crystallographic orientations) within a crystal. An example of a diffusion model fit to a crystal compositional profile is shown in Figure 1c (further profiles are shown in Online Material1 Data A). Our diffusion model best-fits can be grouped into three types:

  1. 39 profiles that show a very good fit with simple modeling and can be simply explained by diffusion with a fixed boundary condition and no crystal growth;

  2. 5 profiles that show a poor fit to the simple model, indicated by an enhanced curvature of the diffusion profile away from the core, relative to the simple model; and

  3. a further 60 profiles that appear to be affected by enhanced curvature (as in 2) and a point of inflection in the profile gradient toward the rim. An example of types 2 and 3 are illustrated in Figure 1c.

All initial modeling accounts for diffusion anisotropy and compositional dependence using the diffusivity data of Dohmen et al. (2007) and Dohmen and Chakraborty (2007a, 2007b). The majority of calculated timescales range from 7 to 45 days (Fig. 4). Outliers (4% of the data set) reach 60–104 days. These longer timescales do not appear to have any relationship with the composition of cores or rims or crystal size. Of particular note is the scatter among the timescales both between crystals and within crystals (Fig. 4), i.e., even accounting for anisotropy and composition dependence. Individual crystals do not consistently correct to a singular diffusion timescale.

Initial discussion

Scatter within the diffusion timescale data set

Scatter between and within crystals could at least partly be due to each individual crystal petrogenesis. It cannot be assumed that the diffusion clock in each crystal was set synchronously and so it is not expected that they should form a timescale “event,” as has often previously been inferred in prior diffusion studies. Since each profile is already corrected for anisotropy, the scatter within crystals is therefore likely to be attributable to sectioning effects, simultaneous crystal growth, and changing boundary condition effects and out-of-plane-diffusion effects (Shea et al. 2015a), as discussed in the following sections.

Sectioning angle

The ideal crystal sections to study for retrieval of accurate diffusion timescales are those that cut through the center of a crystal and also perpendicular to at least one crystal face [and therefore also perpendicular to the diffusion front (Pearce 1984)]. Thin sections cut such that crystals are not sectioned in this way (giving shallow sections) will lead to an apparent lengthening of the observed profiles (Costa and Chakraborty 2004; Costa et al. 2008; Costa and Morgan 2010; Shea et al. 2015a), which is an assumption that is often not explicitly declared or corrected for. This artificial lengthening contributes to an epistemic uncertainty on the timescale introducing a stretch factor to the traverses, most often inducing positive scatter into the timescale data set.

Choosing the narrowest profile within a crystal should be favored to minimize this source of error (Costa et al. 2003, 2008; Costa and Dungan 2005; Costa and Morgan 2010; Shea et al. 2015a). In practice, however, this is not always possible. Crystals picked from tephra and set in grain mounts are often fragmented, and sometimes only one crystal margin is observed, so no relative comparisons can be made. Features such as melt and mineral inclusions, cracks, polishing defects, sub-grain boundaries, or embayments in the crystal considered may further limit options for extracting a profile.

Pearce (1984) calculated the probability of obtaining an ideal cut section of a crystal for the study of zoned crystals. The probability of a section to be within 40% of the center of a crystal that is also within 10° of perpendicular to any one of the three major faces would be around ~20%. This implies a significant proportion of sections are likely to not meet these criteria and be affected by off-center or shallow sectioning. Better representation may be expected, however, when considering larger crystals with a narrower rim thickness, akin to the PdlF olivines in this study. Since sectioning effects are rarely measured and incorporated into diffusion data sets [e.g., Martin et al. (2008) via U-stage measurements] this factor clearly carries some potential as a source of uncertainty.

3D texture

Out-of-plane features that may influence the development of diffusive chemical profiles include proximity to other crystals or vesicles, or low-angle adjacent crystal faces removed during the thin section process where cuts are quite off center and almost graze a face. Crystal morphology may lead to isolated and apparently distinct crystals in thin section that in reality are connected to parts of the same crystal (Welsch et al. 2013). Such crystals might be separately analyzed, and therefore be over-represented in a subsequent data set. Without 3D context, these scenarios represent another source of epistemic error at the scale of single-crystal histories that deserves consideration, which is becoming more widely adopted as a focus for research (e.g., Jerram et al. 2018).

Simultaneous crystal growth and changing boundary conditions with diffusion

Profiles with an enhanced curvature of the diffusion profile away from the crystal core, relative to the simple model, and that also exhibit an inflection point in the profile gradient near the crystal rim (e.g., that in Fig. 1c) are suspected of being affected by crystal growth and/or changing boundary conditions simultaneous with diffusion. These profiles would normally be rejected from further consideration as they do not meet quality control criteria for the degree of match. For instance, the rejection rate would be ~60% for this data set (Online Material1 Data A). These profiles, however, contain valuable information with regards to the environment in which they resided, i.e., a change in ambient conditions. Their rejection, therefore, constitutes a bias, whose impact on both the chronometry and petrographic interpretation is typically ignored. This study attempts to recover this information through dynamic modeling.

Advanced methods

Calculating shallow sectioning effects

The stretching of a compositional profile due to shallow sectioning, as illustrated in Figure 5, can be corrected for using the following equation (Costa and Morgan 2010);


where lm is the measured traverse length, lt is the true traverse length, and the angle θ is the angle between the crystal boundary and the normal to the sectioning plane—the “hade” (Costa and Morgan 2010). We measure this angle using a Universal Stage (U-Stage), the operation of which is described in Kile (2009). The shallower the cut, and hence the larger the value of θ, the more significant the stretch factor on the timescale.

The large size of the olivine grains we used for this study relative to the diffusion scale means that sections could be considerably off center without varying the crystal face intersection angles. Stretch factors were calculated for our data set using the angles measured with the U-stage. At a population level, this serves to highlight how shallow sectioning combined with traditional in-section-plane diffusivities would affect the PdlF diffusion data set overall, and therefore, the possible effects on other data sets. These measurements were repeated to assess uncertainty. Much of the U-stage data set was also re-measured independently by F.K.C. and D.J.M.; the results were found to be reproducible within ±5° by different operators. The stretch factor, stretched timescales, and timescale shift to be applied to timescales modeled by diffusion (that are uncorrected for shallow sectioning) for a range of angles between the minimum and maximum of those measured from this data set can be found in Online Material1 Data A.

Understanding crystal shape using X-ray microtomography

Detailed methodology of the XMT analytical procedure can be found in Online Material1 Data B. Crystals and their textural context were observed in situ in 3D within a rock core using the XMT 3D imaging. We use image segmentation software to segment the olivine crystals from the groundmass within the images produced. To capture, examine in detail, and quantify the shape of the larger (2–3 mm) olivine crystals, several image segmentation methods were tested. Ten slices evenly distributed through the image stack were selected and imaged olivine belonging to large crystals was manually digitized from each slice. The advantage of scrolling above and below the slice (using a z-stack of the image) allowed us to denote olivine as belonging to a large crystal or not, even if the slice intersected only a small area and would otherwise be indistinguishable from groundmass olivine.

“K-means clustering,” “watershed,” and “trainable weka” segmentation modules were applied to these same slices within Image J. K-means clustering partitions a data set into k-groups (clusters). The plugin performs pixel-based segmentation based only on clusters discovered within the image (Wagstaff et al. 2001). The watershed segmentation method considers the input image as a topographic surface, with the brightness of each point representing its height. The entire relief is flooded, and dams or barriers are formed where similar pixels are found, with all similar points forming a catchment basin (Soille and Vincent 1990). The trainable weka segmentation methods combine machine learning algorithms with image processing, assigning areas of interest to “classes” forming the training data set (Arganda-Carreras et al. 2017). The classifier uses a random forest algorithm and decision trees to ultimately classify the selection (Breiman 2001; Hastie et al. 2008).

The brightness and texture of the pixels that make up the olivine are exploited for each segmentation method. Each method was optimized using the 2D slices with manually segmented olivine to choose the best method to be applied to the entire stack. In each case, the result is a binary image that approximates the location of olivine belonging to large crystals. Figure 6 illustrates the good agreement between the percentage area of olivine using all three 3D segmentation methods compared to the percentage area of olivine within the raw 2D slices. We chose to apply the Trainable Weka Segmentation method across the 3D stacks, as this method enabled full automation of each step in the methodology.

Dynamic modeling of crystal growth, changing boundary conditions, and diffusion

To evaluate the potential changes that occur during crystal growth with simultaneous diffusion, we consider two scenarios: (1) that the crystal may grow, and that the outside boundary condition is mobile in space as crystallization proceeds, and (2) that the equilibrium composition of the olivine is changing as crystallization proceeds, due to melt evolution.

In these cases, the growth of a crystal and the change in equilibrium value are driven by changes in magma temperature and the degree of crystallization that this affords. Because diffusivity is a strong function of temperature, we must also consider: (1) that the diffusivity will be controlled by temperature, which will decrease during cooling and crystallization, and (2) the thermal effects on oxygen fugacity.

To constrain this series of processes, we chose to model diffusion over the crystallization path of the system as a function of temperature, to obtain the equilibrium olivine values at any given temperature, which can be applied as a boundary condition. By specifying a cooling rate applied to the sample, this can be expected to control the evolution of the external melt composition with time and control the diffusivity as a function of time. What the diffusion modeling requires is a precise starting melt composition, which can be modeled using software such as Petrolog3 (Danyushevsky and Plechov 2011) or MELTS (Ghiorso and Sack 1995; Asimow and Ghiorso 1998). We used Petrolog3 to determine the initial melt composition prior to olivine rim growth and groundmass formation. We extracted variable amounts of Fo84 olivine from the whole-rock composition and assessed olivine-melt equilibrium. We then project crystallization of this melt considering the phases present in the observed assemblages (with the exception of orthopyroxene, which is present in low abundances in the Petrolog3 models). A detailed methodology is described in Online Material1 Data B.

Thermodynamic modeling using Petrolog3.

Olivine in the models typically crystallizes through the interval Fo84–79, but more evolved olivine is not present; it is replaced in the models by orthopyroxene, which crystallizes at the level of a few percent. This is distinct from what we see in the specimens, where no orthopyroxene is observed, and olivine rims exhibit compositions down to Fo65 or less (indicated by brighter BSE pixels at the very rim, below the spatial resolution of the EPMA). We suspect this represents a limitation in the modeling vs. nature; either the nucleation of orthopyroxene is somehow suppressed in the natural samples or that the crystals are so small as to be visually unidentifiable from olivine within the fine-grained groundmass. Regardless, orthopyroxene does not form overgrowths on olivine in RU0701, and as such there is no evidence to support its inclusion across the petrogenetic interval of interest. Therefore, we excluded orthopyroxene from the Petrolog3 modeling, restricting the assemblage to olivine, plagioclase, clinopyroxene, and oxide. Olivine crystals in natural samples were able to exchange via diffusion with an adjacent, more-evolved melt that equilibrate with Fo65, since fractionation to a cumulate has not occurred.

Figure 7 shows how crystallization behavior changes with temperature in the model. This curve is parameterized for use in a diffusion model—each stretch of cotectic is fitted with a third-order polynomial with R2 > 0.99, which is used within that temperature range to determine the equilibrium olivine composition for any given temperature. Note that at temperatures below 1117 °C, olivine does not crystallize within the model, and the values reflect the equilibrium due to exchange with an evolved ambient melt and a Kd for olivine of 0.303.

Dynamic diffusion modeling.

The modeling was conducted using a simple, iterative one-dimensional finite difference model, whose logic follows that shown in Figure 8. In this manner, the model allows for tracking the liquidus and an evolving melt, allowing the temperature to control both the composition, the diffusivity, and the absolute oxygen fugacity (via a buffer), much as it would during the cooling of the magma.

Advanced results

Sectioning effects

Our analysis of sectioning effects shows a large frequency of sections cut within angles 10° of vertical and a smaller frequency of the very shallow angles (where the crystal edge is sectioned up to 30° to the vertical; Fig. 9). These angles correspond to timescale shortening between 0–2.5 and 33.3% relative. The average angle measured from vertical (the “sectioning” angle) for the PdlF data set is 8.9° and would require an average shift of ~2.5% applied to the timescale data set (0.011 log units), which is low, relative to other sources of uncertainty. A full list of measured angles and the corresponding shift on individual timescales is included in Online Material1 Data A.

Quantification of 3D crystal shapes

3D renders of four olivine crystals or crystal aggregates are shown in Figure 10. A further 25 crystals are provided in Online Material1 Data A. Assessment of these images reveals that only 25% of the crystals could be considered euhedral; 68% are subhedral, and 7% are anhedral. Furthermore, only 25% are individual crystals and the remainder are present as crystal agglomerates. These agglomerates are generally crystallographically aligned (see Fig. 10e). Nearly all the crystals are polyhedral, apart from two that are tabular shaped; 19% are rounded nodules, and 13% of the agglomerates show parallel organization or hierarchy as described by (Welsch et al. 2013).

Crystal growth and changing boundary condition incorporated into diffusion timescales

Using the initial melt composition from Petrolog3 (see Fig. 7), we ran several diffusion models (as described in Fig. 8) to assess the effects of changing boundary conditions and crystal growth on diffusion timescales. The results of each are described in the following sections.

Modeling sequence 1: Changing boundary condition (no growth).

Modeling was initially conducted for a starting temperature of 1126 °C, at NNO+1, for various cooling rates. This generates a sequence of diffusion curves with geometries as shown in Figure 11. The profiles follow proportionally similar trajectories through temperature and time, are self-similar and display typical proportionality for diffusion profiles, i.e., where cooling rates quadruple, the diffusion width decreases by a factor of two.

Modeling sequence 2: Changing boundary condition with growth.

Modeling was also conducted for a starting temperature of 1126 °C, at NNO+1, for a single cooling rate of −1 °C/h, combined with various growth rates. This generates a sequence of diffusion curves with geometries as shown in Figure 12; the addition of growth stretches the profiles, as would be expected. As growth is considered here to be linear in time, and diffusion is reducing exponentially due to linear cooling in time, the profiles are decoupled and are not self-similar. The highest growth rates are in effect recording a dominance of the liquidus surface, which is not completely smooth, and which has been subjected to a variable degree of diffusional smoothing with distance.

Modeling sequence 3: Fitting a specific case.

Due to the dynamic interplay of parameters, an initial profile was considered for detailed investigation using these diffusion and growth models. The fitting was performed iteratively by changing the relative magnitude of cooling rates and growth rates to gradually improve the fit for a particular traverse. To allow for an appreciation of both the similarity and difference between the curves in Figures 11 and 12, and real data, we superpose a profile onto these curves (labeled “measured”).

These figures show that natural data are significantly divergent from both cases studied, though growth allows us to fit the broader rim. A particular feature of the natural curve is that it has a gradual concave curvature, while the diffusion-only curves are convex. Such concavity is consistent with increasing growth rates; this can be seen on Figure 12, as the best fit seems consistent with high growth rates at the rim, but with lower growth rates near the core. Crystals, therefore, seem to have experienced a late-stage, accelerating rim growth. This was considered by applying a temperature-dependent growth across a temperature range such that:


where G is the growth rate, G0 the half-growth rate at Gmin, T0 is the start temperature and Tmin is a nominal finish temperature. This produces a gradually increasing growth rate following a parabolic curve. By adjusting the parameters of cooling rate and the parameter G0, above, we can attempt to emulate the curvature on the natural profile.

Model solutions can attain quite good agreement with the natural data using this method (Fig. 13). Slight deviations near the core may suggest a degree of two-step history, with some degree of diffusion prior to surface emplacement, although to explicitly fit this is hard to justify. Scenarios involving the departure of an individual crystal from the single-cooling trend for a period of time are not difficult to imagine in nature, yet we adhere to the cooling and growth rates expected from the system scale observations summarized by the model (Fig. 7).

Extension to the broader Piton de la Fournaise data set.

Using the dynamic methods, we re-modeled 13 profiles from the RU0701 data set. We selected profiles that exhibited a poor fit when initially modeled at a single temperature of 1126 °C and fixed boundary condition.

By incorporating a cooling rate and a growth rate guided by the expected behavior of the cooling magma (Fig. 7) the model produced results that are in good agreement with the natural profiles. Six of these profiles from two crystals (RU0701_2_G and RU0701_2_I) are shown in Figure 14. The remaining six re-modeled profiles are included in Online Material1 Data A. All of the re-modeled profiles are consistent with an average cooling rate of −0.5 °C/h (±0.1 °C) and an average growth rate of ~1.9 × 10–11 ms–1, over the cooling range 1127–967 °C. These models translate to new diffusion timescales of, on average, 13 days (with a minimum timescale of 10 days and a maximum timescale of 17 days).

The textural evidence, profile shapes at the crystal rims, and short diffusion timescales indicate that the diffusion rims are late-stage and most likely related to lava flow cooling. As an external check on the viability of our re-modeled timescales, we would expect the lava flow cooling time and the re-modeled diffusion timescales to be similar. We can use the following equation to calculate the cooling time of the lava flow:


where t is the time in seconds, A is the thickness of the lava flow in meters, and K is the thermal diffusivity of a basalt (Hon et al. 1994). This simple lava flow cooling rate calculation (with a lava flow thickness of 2–3 m) gives cooling rates of 9–19 days. The independent line of cooling rate calculation agrees with the remodeled diffusion timescales using the models that incorporate cooling rate (and thereby changing boundary condition) and crystal growth rate, and so lends confidence to our approach.

Advanced discussion

It is reasonable to consider that the diffusion “clock” in each crystal can be set at different times yet stop at effectively the same time (on eruption). Scatter within crystals due to natural variation, that if understood and placed into broader petrological context, can lead to detailed understanding of magmatic and/or volcanic phenomena (e.g., Pankhurst et al. 2018b). Diffusion timescale scatter between crystals could, however, also occur even if diffusion clocks are set simultaneously if complicating factors such as proximity to other crystals, exposure to different melts, mineral grains, or vesicles blocking diffusion have persisted for a significant fraction of time and if simplistic modeling is applied. Scatter within individual crystals could reflect limitations on the quality of modeling and/or be due to real complications such as crystal growth, changing boundary conditions, sectioning effects, and 3D crystal geometry.

Shallow sectioning effects

There is limited control for ideal alignment of crystals during the preparation of grain mounts, and effectively no control in the case of thin sections. But in contrast to the expected bias toward shallow angles (Pearce 1984), our data set largely consists of boundaries with relatively small deviations (<10°) from being ideally sectioned. As a consequence, accounting for the measured sectioning angle has a relatively limited effect of reducing the timescales from mostly 7–45 days (outliers between 60–104 days) to 6–44 days (outliers between 62–82 days; Fig. 15).

Since the combined uncertainty of binary Fe-Mg inter-diffusion in olivine due to temperature, activation energy, and D0 is 0.34–0.38 log units, an average 0.011 log shift toward longer timescales cannot be considered critical at a population level, for olivine grains of this large size relative to their narrow diffusion width.

Overall, the relatively low occurrence of significantly shallow sectioning angles within the RU0701 data set is likely due to the unconscious selection of crystal slices and faces analyzed combined with the large relative difference in crystal size to profile length. Even when trying to remove bias, unworkable crystal slices do not make the initial selection. Sectioning angle uncertainty is likely to become more important in small data sets and in circumstances where context is lacking or where only one crystal margin can be considered.

Understanding crystal shape

Shea et al. (2015a) conducted a comparison between 1D, 2D, and 3D diffusion modeling of timescales from slices based on theoretical polyhedral, spherical, and cuboid shapes. They found that timescale distributions vary depending on the crystal faces present and the angles between the faces. These authors also showed that the difference between measured/calculated to real timescales can be especially variable where diffusion fronts from adjacent crystal faces converge, i.e., where profiles are extracted nearer corners. This results in a longer apparent profile and therefore an apparent longer timescale when modeled. These zones are a common feature of polyhedral crystal slices (a common natural crystal shape) where faces meet at <180°.

We considered the Piton crystal form, as the rounding of crystal faces will enhance the effect of merging diffusion fronts and the arrangement of crystals. We also considered if grains were single crystals or present as aggregates, and how these came together, whether by syneusis (Schwindinger and Anderson 1989) or paired nucleation. Diffusion timescales may be affected by a “pre-history” of individual parts, a change in boundary conditions that coincide with their coming together.

An observation made possible here due to the 3D data is the absence of an (001) face; the theoretical olivine shape used in the study of Shea et al. (2015a) does have an (001) face. Furthermore, the 3D data demonstrate that the olivine crystals do not conform to an “ideal” solitary crystal morphology. An initial comparison of a slice through the polyhedral olivine shape used by Shea et al. (2015a) and a slice through a 3D rendered natural crystal at the same orientation is sufficient to illustrate that the natural crystal is not perfectly euhedral and that additional faces are present within crystals (Fig. 16). In this case, a timescale calculation is unlikely to suffer much uncertainty due to the relatively narrow margins and wide areas from which to choose a profile away from corners. But without 3D data, this view would only be intersected by chance, and in the cases of smaller crystals with broader margins even subtle departures from an ideal shape will magnify the geometric uncertainties when studying natural crystals. Further work, and more comprehensive natural 3D data, is required to better link the value gained from using generalized synthetic crystal shapes with applications across various settings.

The XMT data set shows that natural crystals have added morphological complexity that may interact with diffused zones to a greater extent than has been considered in previous studies where only idealized morphologies have been investigated. It is clear that the effects of euhedrality and arrangement of the crystal agglomerates need to be more fully quantified to understand how crystal shape uncertainties may affect calculated diffusion times-cales. Uncertainty on calculated diffusion timescales would likely increase, yet placing “hard” numbers on this effect is difficult to establish since we used a polychromatic beam prior to method developments that attempt to address the inherent problems in deriving quantitative attenuation data using polychromatic beams (Pankhurst et al. 2014, 2018b, 2018c).

Fully quantitative 3D composition data throughout olivine crystals, which at this stage of research is still only possible using a monochromatic X-ray beams, is needed to address these questions. Current work strives to develop a 3D method that uses polychromatic beams, is rapid, chemically quantitative, and able to resolve high-frequency gradients such as in crystal margins, with the same confidence as 2D imaging. At present, the 3D data provides the important observations that this natural crystal population is strongly affected by crystal rounding and is mostly present as crystal agglomerates.

Simultaneous crystal growth and changing boundary conditions

Each of the 13 dynamically modeled profiles were remodeled with a global cooling rate, and as the magma cools, the boundary condition continuously shifts, with simultaneous crystal growth. The cooling rate used for each profile (–0.5 ± 0.1 °C/h) is consistent with diffusion related to quenching and lava flow cooling and inferred from textural observations (an average peak growth rate of ~1.9 × 10–11 ms–1 was applied to each profile). Most profiles (11 of 13 profiles) were re-modeled using an accelerating growth rate, rather than a linear growth rate, as the growth rate seems also to be a function of temperature and has to be increasing with continued cooling in most cases, based on the profile shapes.

Although the profiles show an apparent increased diffusion length due to crystal growth and changing boundary conditions effects, our advanced models accommodate how diffusion speed slows as the lava flow cools, and so the diffusion timescales do not become significantly shorter nor is the scatter of the data set significantly reduced. As a result of reducing epistemic uncertainty, there is an attendant increase in confidence that the scatter reflects real variation in crystal histories. This is corroborated by the observation of different final rim compositions (between Fo72 and Fo60): even though the crystals reflect a common cooling environment, the details of their final experiences were controlled locally by the availability of melt and other considerations, including pore space. This appears to be particularly important in the case of samples from lava flows where post-eruption cooling must be considered in petrogenetic interpretations, which could, in turn, free models of deeper magmatic plumbing system structure and dynamics from this frequently observed complexity. Cheng et al. (2020) suggest that it can be possible to create simple diffusion curves (as seen within our samples) as a result of changing boundary conditions due to flow dynamics if there is sufficient residence time for complex zoning patterns to evolve into simple zoning—we do not consider this to be the case for these PdlF olivines.


Uncertainty of individual diffusion timescales can be considerably reduced using a combination of three-dimensional observation and dynamic forward modeling that includes an assessment of growth and changing boundary conditions. The quantification of the sources of uncertainty upon the RU0701 data set allows us to better understand how they, and by extension, their corrections, affect timescale distributions and interpretation. The results have implications for the design of crystal diffusion studies, which we recommend should explicitly consider sectioning angle, 3D geometries, crystal growth, or changing boundary condition corrections on reported timescales. A re-appraisal of modeling assumptions, as demonstrated here, may significantly improve both the quality and quantity of data recovered from a suite of crystals. It follows that increased data set size per investment in resources, as well as increased confidence in timescale outputs, will have value as kinetic and kinematic views of magmatic processes improve toward full physical history reconstruction (e.g., Pankhurst et al. 2018b).

Our case study demonstrates that agreement can be found between the timescales of cooling at the scale of a lava flow and the range exhibited by the individual crystals it contains. This result would not have been possible without the need to consider dynamic modeling as a solution to poorly fitting initial profiles to simplistic diffusion models. The system-scale information required, such as intensive parameters and a model of melt evolution with cooling, in turn, provide a framework to better contextualize the variations observed in final rim composition and timescales. The consistently high-quality fits between raw data and models using a common cooling history and growth rate lead us to frame new questions regarding the variation of rim compositions and other details, which will be addressed in a subsequent contribution. It also provides a demonstration that the cooling regime within a lava flow can be obtained on the crystal scale, which likely holds value in understanding more complex flow settings.

The value in improving the quality of diffusion timescale data sets is expected to provide a better return on data-gathering investment, including observations at the population scale that are based on larger and more complete data sets and increase the detail and therefore utility of interpretations.


F.K.C. and D.J.M. acknowledge NERC Studentship number 1367441 for supporting this work. M.J.P. acknowledges an AXA Research Fund Fellowship that supported part of this work. Sample collection was made possible by the National Geographic Society (NGS 8330-07). XMT scanning was funded through UK-NERC grants NE/M013561/1 and NE/N018575/1.


We thank the reviewers and editor for their constructive comments and the handling of the manuscript. We thank T. Shea for the helpful discussions and use of the theoretical olivine shape MATLAB code. J.W. Williams and R. Walshaw are thanked for their assistance with sample preparation and the analytical facilities at the University of Leeds. R. Clark and A. Booth are thanked for providing valuable writing time to enable the preparation of this work for publication. L. Courtois and S. Nonni are thanked for assistance with XMT scanning.


Deposit item AM-21-37296, Online Material. Deposit items are free to all readers and found on the MSA website, via the specific issue's Table of Contents (go to http://www.minsocam.org/MSA/AmMin/TOC/2021/Mar2021_data/Mar2021_data.html).
This is an open-access article distributed under the terms of the Creative Commons Attribution CC-BY 4.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.Open access: Article available to all readers online. This article is CC-BY.