# Geology

## Abstract

A random search modeling approach is introduced to retrieve the thermal histories of magmatic bodies from chemical and isotopic analyses of olivine crystals. The sample used for demonstration is an olivine phenocryst from Kilauea Iki lava lake, Hawaii. Random time-temperature (*t*-*T*) paths were used to model diffusion accompanying crystal growth. Although a large number of *t*-*T* paths yielded acceptable fits to the Mg/Fe chemical profiles measured in olivine, only a small fraction resulted in acceptable fits to the Mg-Fe isotopic profiles. The best-fit cooling path was close to the independently inferred cooling history of this sample. Use of isotopic profiles in minerals coupled with a random search method can be a powerful tool for recovering thermal histories of rocks with little geologic context.

## INTRODUCTION

Chemical zoning in olivine crystals records a wealth of information on the cooling and crystallization history of magmatic bodies. However, interpreting the nature of chemical zoning is a challenge because chemical diffusion and crystal growth can produce the same chemical zoning patterns. Misinterpreting zoning profiles produced by crystal growth as arising from diffusion would lead to calculated diffusion time scales that are meaningless.

In most natural systems, diffusion and crystal growth are expected to act concurrently, as diffusive fluxes arise from changes in magma composition induced by crystallization. Studies have shown that negatively coupled Mg-Fe isotopic profiles are associated with chemical diffusion as these elements interdiffuse, and the light isotopes always diffuse faster than their heavier counterparts (Dauphas et al., 2010; Teng et al., 2011; Sio et al., 2013; Oeser et al., 2015). If chemical zoning is controlled by crystal growth, little Mg and Fe isotopic fractionation is expected because the equilibrium isotopic fractionations between olivine and silicate melt are small at temperatures relevant to magmatic systems (Teng et al., 2007; Dauphas et al., 2014). Modeling of the isotopic profiles is thus a powerful approach for discerning the contribution of diffusion versus growth in generating chemical zoning.

To fully harvest the potential of chemically zoned minerals to record the thermal histories of magmatic bodies, we developed a random search modeling approach that calculates the time-temperature (*t*-*T*) histories consistent with the measured chemical and isotopic profiles. This approach has been tested successfully on a sample from the Kilauea Iki lava lake, and it can be extended to a wider variety of geological settings.

## SAMPLE

The olivine phenocryst studied by Sio et al. (2013; see also Fig. DR1 in the GSA Data Repository^{1}) was used to test the modeling approach because its cooling history is well constrained. This sample (KI81–5-254.5) came from Kilauea Iki lava lake, Hawaii, which was formed in A.D. 1959 by inflow of magma into a preexisting crater pit. The eruptive magma temperature was ∼1230 °C based on the liquidus temperature and the core composition of olivine phenocrysts (Helz et al., 2014). Upon eruption, the lava cooled to ∼1190 °C, from which point it cooled more slowly in the lava lake. The thermal history for this sample was perturbed by nearby core drilling activities in 1981, resulting in rapid cooling from ∼1130 °C to ∼1070 °C in the 3 wk prior to retrieval of the sample by the U.S. Geological Survey (USGS; Helz and Wright, 1983; Sio et al., 2013). The olivine phenocryst studied has a rim of overgrowth due to this late-stage rapid cooling episode.

Sio et al. (2013) demonstrated the presence of large and negatively coupled Mg-Fe isotopic fractionations in this olivine phenocryst. As a first pass, they used a simplified no-growth model and discarded the late-stage overgrowth to estimate the isotope effect of diffusion for Mg and Fe in olivine. In the present paper, crystal growth is incorporated to enable rim-to-rim modeling of the chemical and isotopic profiles.

## METHODS

### Monte Carlo Inversion

Our random search approach is similar to the Monte Carlo inversion method used in geophysics (Press, 1968) and the inverse-forward modeling technique used in thermochronology (e.g., Ketcham, 2005). Geologic constraints are represented as Monte Carlo boxes in a time-temperature (*t*-*T*) space. The *t*-*T* trajectories are then created by connecting points that are chosen randomly within these boxes, as explained in more detail in the following.

Although the thermal history of our sample is well constrained by petrologic studies, we used modeling parameters calculated from the MELTS software (melts.ofm-research.org; Ghiorso and Sack, 1995) in order to demonstrate the versatility of the random search approach. Using a best estimate of the starting composition of the magma (Table DR1), the liquidus temperature calculated using MELTS is 1230 °C, and this is the assumed initial temperature of the magma. After assigning a ±10 °C error bar and defining time zero to be the formation of the lava lake, the first constraint in the (*t*-*T*) space is (0, 1230 ± 10). The composition of olivine in equilibrium with the melt at any temperature is calculated using MELTS. The quench temperature is taken to be 1085 °C, corresponding to the measured rim composition of the olivine (Fo_{68}). Hence, our second constraint is placed at (*t*_{q}, 1085 ± 10) in the (*t*-*T*) space, where *t*_{q} is time of quench (= 21.4 yr).

With these constraints, trajectories are created such that the first point (*t*_{1}, *T*_{1}) is at (0, 1230 ± 10). Assuming monotonic cooling of the lava lake, the second point (*t*_{2}, *T*_{2}) is chosen randomly with equal probability in the space occupied between (*t*_{1}, *T*_{1}) and the end point (*t*_{q}, 1085 ± 10). The third point is again chosen between (*t*_{2}, *T*_{2}) and the end point (*t*_{q}, 1085 ± 10) and so forth. Each *t*-*T* path is defined by 10 points, which is more than sufficient to capture the diversity of possible cooling paths (Fig. 1). Sixty (60) random *t*-*T* trajectories are generated. Forward modeling is done for every trajectory. The calculated and measured chemical profiles are then compared, and isotopic profiles are calculated only when there is a good match to the chemical profile.

### Forward Modeling of Diffusion Accompanying Crystal Growth

The calculations are done using the nondimensional equations given in the GSA Data Repository. The nondimensional diffusion equations are solved using the explicit finite difference method, which has been validated against available analytical solutions for specific initial and boundary conditions (Figs. DR2 and DR3).

The equation governing diffusive exchange in spherical coordinates is (Crank, 1956):where *C* is defined to be the molar fraction of Fe in the octahedral sites in olivine: *C* = Fe/(Fe + Mg), *R* is radius, and *D* is the diffusion coefficient. As Mg-Fe interdiffusion in olivine is anisotropic, we take into account the crystallographic orientation of the olivine when calculating *D* (see the Data Repository). However, this is still an imperfect correction, because Equation 1 is used to describe isotropic diffusion in a sphere (for a rigorous treatment of diffusion anisotropy, see Watson et al., 2010). The effect of imperfectly correcting for diffusion anisotropy can be evaluated by modeling the crystal as having a slab or cylindrical geometry. This in essence represents the scenario where diffusion in some directions is infinitely smaller than others. Figure DR4 shows that incompletely correcting for diffusion anisotropy will likely have an effect of less than 20% on the magnitudes of isotopic fractionation.

Profiles are calculated independently for each isotope (^{54}Fe, ^{56}Fe, ^{24}Mg, and ^{26}Mg). We assume that the composition-dependent diffusivity determined in olivine diffusion experiments represents that relevant to ^{56}Fe and ^{26}Mg interdiffusion. For the lighter isotopes of Fe and Mg, the interdiffusion coefficient is multiplied by a factor equal to the inverse ratio of the masses raised to a β exponent (i.e., the isotope effect for diffusion; Schoen, 1958; Tsuchiyama et al., 1994; Richter et al., 1999): *D*_{24Mg} = *D* × (26/24)^{βMg}, and *D*_{54Fe} = *D* × (56/54)^{βFe}, with β_{Mg} = 0.16, and β_{Fe} = 0.27 (Sio et al., 2013; for comparison, Oeser et al. [2015] reported estimates of 0.055 < β_{Mg} < 0.12 and 0.075 < β_{Fe} < 0.30). The isotopic profiles are calculated as δ^{26}Mg = (*C*_{26Mg}/*C*_{24Mg} – 1) × 1000 and δ^{56}Fe = (*C*_{56Fe}/*C*_{54Fe} – 1) × 1000.

To solve Equation 1, an initial and two boundary conditions are required. The initial condition is the initial Fe concentration (*C*) of the olivine phenocryst, which is the olivine composition at the liquidus calculated by MELTS: *C*(*R*, 0) = 0.145 or Fo_{85.5}, where Fo# = 100 × (1 – *C*). The first boundary condition is given by the symmetry of the phenocryst, which imposes a no-flux boundary condition at the center of the crystal: ∂*C*/∂*R* = 0 at *R* = 0 for *t* ≥ 0. The second boundary condition represents how the olivine compositions at the crystal-melt interface change as a function of temperature. While this boundary condition has been determined empirically (Helz and Thornber, 1987; Helz et al., 2014), it was calculated using MELTS in this study, as we wished to test the proposed modeling approach by assuming that minimal information is known about the sample. The MELTS outputs are given in terms of Fe/(Fe + Mg) for the olivine; parameterizing this ratio as a function of temperature gives *C*(*T*) = (5.3413 × 10^{–6}) × *T*^{2} – (1.3415 × 10^{–2}) × *T* + 8.566, where *T* is in Celsius. This boundary condition defined by MELTS is close to that independently determined for Kilauea Iki lava lake samples (Fig. DR5), suggesting that uncertainties in the starting magma composition and thermodynamic treatment of Mg-Fe partitioning between olivine and magma in MELTS do not influence the boundary condition significantly.

The position of the boundary condition at the moving interface depends on the crystal growth rate. Olivine growth rate in basaltic melts is poorly constrained. Experimental studies for olivine growth with small degrees of undercooling have only been done for short durations (Jambon et al., 1992). Cashman (1993) found a positive relationship between plagioclase growth rate and cooling rate in samples from Makaopuhi lava lake, Hawaii (Fig. DR6). We assume that this roughly linear relationship also exists for olivine and thus parameterize the olivine growth rate as follows: υ = ∂*R*/∂*t* = –γ × *dT*/*dt*, where γ is an input and has units of µm/°C. While this functional form for growth rate is highly uncertain (see the Data Repository), it nonetheless provides a first-order approximation that can be improved upon once the parameters governing olivine crystal growth in natural settings are better understood. The effect of crystal growth on isotopic fractionation is shown in Figure DR7.

The moving boundary position *R*_{B}(*t*) is thus:

As the olivine grows, we assume no equilibrium (or surface reaction–controlled kinetic) isotopic fractionation of both Mg and Fe isotopes between olivine and melt. For Mg, this is inferred from the lack of Mg isotopic fractionation during magmatic differentiation (Teng et al., 2007). For Fe, synchrotron measurements of equilibrium isotopic fractionation between olivine and basaltic glass suggest that δ^{56}Fe_{melt} – δ^{56}Fe_{olivine} ≈ 0.03‰ ± 0.03‰ at 1130 °C (Dauphas et al., 2014; assuming Fe^{3+}/Fe_{total} = 0.14 in the melt).

To assess the effects of uncertainties in diffusivity, crystal growth rates, and β values, we run models varying these parameters. For diffusivity, we do this by applying a multiplication factor: *D* = κ× *D*_{traverse}, where *D* is the diffusivity used in the models, κ is the multiplication factor, and *D*_{traverse} is the diffusivity calculated along the crystallographic direction of the analyzed traverse. Values of 1, 2, or 3 are assigned to κ, because higher diffusivities than those calculated by the parameterization of Dohmen and Chakraborty (2007) have been reported (e.g., Misener, 1974; Jurewicz and Watson, 1988).

For growth rates, we assign values of 3, 4, or 5 µm/°C to γ. The integrated growth in the model is calculated from Δ*R* = γ (*T*_{0} – *T*_{q}). The values chosen for γ yield Δ*R* of ∼400–700 µm, which seems reasonable based on the size of the isotopically flat region near the rim of the crystal, which probably represents overgrowth during the late-stage rapid cooling episode (Fig. 2; Fig. DR1).

Last, isotopic fractionations scale linearly with β values so that the modeled Mg and Fe isotopic fractionations can easily be calculated for values other than those adopted by Sio et al. (2013; β_{Mg} = 0.16 and β_{Fe} = 0.27). Results are also shown for β_{Mg} = 0.12 and β_{Fe} = 0.23 (keeping β_{Fe}/β_{Mg} ∼2; Sio et al., 2013; Oeser et al., 2015). Parameters used and modeling results are given in Table DR2.

### Evaluation of Goodness of Fit

The global best fit is found by searching for the smallest reduced chi-square, , across all models. After identifying the global best fit, an unknown (geologic) source of error, σ_{unk}, is added to make its equal to 1. The same σ_{unk} is then applied to all models:where σ_{i} is the associated uncertainty on the isotopic analysis at position *i* across a profile, *v* = *N* – *n* – 1 is the number of degrees of freedom, *N* is the number of data points, and *n* is the number of parameters. In our case, given the correlation between δ^{26}Mg and δ^{56}Fe, the number of data points is the number of positions where isotopic data are taken (*N* = 8). The number of parameters is the number of nodes in the cooling history. Although eight nodes are used in each trajectory, there are typically five unused inflection points (i.e., points that do not influence the *t*-*T* path because they merge with the quench point). For three effective inflection points, we have *n* = 6, describing randomness in both time and temperature. Diffusivity, growth rates, and β values are not considered parameters because these values are set (i.e., they are not fit parameters). Acceptable isotopic fits are those that have ≤ 3.84, denoting a one-sided significance level of 0.05 for 1 degree of freedom.

## DISCUSSION

We consider that any model able to reproduce the measured core composition (Fo_{83}) of the olivine phenocryst to within ±0.5 in Fo# is satisfactory. Isotopic profile calculations are performed only when satisfactory chemical profiles have been obtained. The global best fit is found using the parameters γ = 5 µm/°C, *D* = 2 × *D*_{traverse}, β_{Mg} = 0.16, and β_{Fe} = 0.27 and is highlighted in green in Figure 1. The acceptable fits, defined as the ones with ≤ 3.84, are colored in yellow, and they cluster closely around the best-fit *t*-*T* trajectory (Fig. 1). Acceptable fits using other sets of parameters also resulted in *t*-*T* trajectories that are identical to the global best fit, demonstrating that this modeling approach is insensitive to reasonable uncertainties in *D*, γ, and β values (Figs. DR8–DR12).

Many of the *t*-*T* trajectories that reproduce the chemical profile fail to reproduce the Mg and Fe isotopic profiles. The *t*-*T* trajectories associated with the acceptable fits are most similar to the independently inferred *t*-*T* history of the lava lake (in open circles in Fig. 1). Particularly, the late-stage rapid cooling experienced by the sample is well captured by this random search modeling approach. This demonstrates that the isotopic data provide constraints on the cooling histories of magmatic bodies that are not provided by chemical profiles alone. Chemical profiles alone are nondiscriminating in that they cannot detect the late-stage overgrowth that extends the zoning profile of the interior. Isotopes, on the other hand, are not very fractionated in this overgrowth, indicating that the sample experienced late-stage rapid cooling.

The *t*-*T* trajectories can be scaled with time, and acceptable fits can still be found using different sets of parameters (*D*, γ, and β values). Figures DR13 and DR14 show the same *t*-*T* trajectories but scaled to quench times (*t*_{q}) of 12 yr and 60 yr. Using β_{Mg} = 0.16 and β_{Fe} = 0.27, the parameter sets γ = 5 μm/°C, *D* = 3 × *D*_{traverse} and γ = 3 μm/°C, *D* = *D*_{traverse} can be used to fit profiles with *t*_{q} = 12 yr and *t*_{q} = 60 yr, respectively (Figs. DR13 and DR14). These are the shortest and longest cooling time scales allowable by the ranges of *D* and γ explored in this study. Therefore, while this random search modeling approach is still susceptible to uncertainties in *D*, γ, and β values in retrieving absolute cooling time scales, the morphology of the retrieved cooling path is robust. That is, even though we vary *D*, γ, and β values, and *t*_{q} in our calculations, we consistently obtain *t*-*T* trajectories that describe a late-stage episode of fast crystal growth accompanied by little diffusion.

Recently, Oeser et al., (2015) used β_{Mg} from 0.055 to 0.12 and β_{Fe} from 0.075 to 0.3 to fit isotopic profiles measured in olivines from Massif Central, France, and the Canary Islands, Spain. This wide range of β values may arise from uncertainties in boundary conditions for their diffusion models. In addition, β values could possibly depend on crystal composition, crystallographic orientation, and temperature (Van Orman and Krawczynski, 2015). Therefore, much complexity may be hidden in the effective β values used in this study. Laboratory experiments to study the mechanisms controlling β values are necessary to apply this chemical-isotopic modeling approach to samples with less geological contexts.

## CONCLUSIONS

The random search modeling approach proposed here can be used to retrieve thermal histories of geological samples. This approach is superior to using the chemical profiles alone because isotopes are sensitive to diffusive processes, allowing chemical zoning that arises from diffusion to be distinguished from crystal growth in an evolving magma. There are still uncertainties in the values of Mg-Fe interdiffusivity in olivine, the isotope effect for diffusion, and the growth rate of olivine that limit the precision of reconstructed cooling time scales. Nonetheless, we demonstrate that by combining petrographic, chemical, and isotopic studies of zoned minerals, much tighter constraints on thermal evolutions can be achieved. The modeling approach will be particularly useful in volcanology for understanding magmatic processes prior to and during eruptions.

## Acknowledgments

Frank Richter, Jonathan Gagné, and Neil Bennett are thanked for insightful discussions. Helpful reviews by James M. Watkins, E. Bruce Watson, and Megan E. Newcombe greatly improved this manuscript. This work was supported by a Carnegie Postdoctoral Fellowship and a NASA Earth and Space Science Fellowship to Sio, and the National Science Foundation (EAR144495, EAR150259) and NASA (NNX14AK09G, OJ-30381–0036A, NNX15AJ25G) grants to Dauphas.

## Footnotes

↵

^{1}GSA Data Repository item 2017017, numerical formulations, discussions on the functional form of crystal growth rate, and additional plots showing how uncertainties in modeling parameters affect the calculated time-temperature histories, is available online at www.geosociety.org/pubs/ft2017.htm, or on request from editing{at}geosociety.org.

- Received 2 May 2016.
- Revision received 17 October 2016.
- Accepted 17 October 2016.

- © 2016 Geological Society of America