Phosphorus Coupling Obfuscates Lithium Geospeedometry in Olivine

Lithium zoning in zircon and olivine provides a powerful record of crystal histories and magmatic processes in volcanic systems. Characterizing Li behavior in olivine is important because it is one of the few elements in basaltic systems that track short-duration magmatic processes (e.g., less than a few days). However, the potential for trace element coupling and the competing effects of crystal growth and subsequent diffusion obfuscate interpretations of Li zoning patterns and their use for geospeedometry. Here, we use diverse analytical techniques (EPMA, LA-ICPMS, nanoSIMS) to untangle records of growth and diffusion in carefully oriented olivine crystals from Kīlauea (Hawai‘i). Li, P, and Al are targeted because they show correlations despite contrasting behaviors during growth and diffusion. Lithium zoning exhibits two styles: (1) non-coupled, wherein Li (1–3 ppm) shows diffuse zoning over 10s–100s of μm between the crystal core and rim with low (10s of ppm) and homogeneous P, and (2) coupled, wherein sub-ppm Li enrichments <40 μm wide are correlated with P-rich zones (100s of ppm) that preserve a blueprint of crystal growth. Non-coupled Li zoning modeled for diffusive re-equilibration yields hours–days timescales, reflecting magma-mixing events that primed the system for eruption. Coupled Li peaks consistently occur where analyses traverse the dendritic framework of P zoning that reflects rapid crystal growth. This correlation results from Li+ ions partially satisfying charge balancing requirements of P5+ ions in the olivine lattice. Aluminum balances most of the P budget in olivine but has non-systematic correlations with Li peaks. Transects can exhibit both non-coupled and coupled Li behaviors, preserving histories of both crystal growth and diffusive re-equilibration. Modeling a near-rim Li enrichment peak with the Li diffusion coefficient yields a timescale of <3 min, too short to represent the time elapsed following crystallization. If the near-rim Li and P peaks are modeled using the diffusion coefficient for P, congruent timescales of 5–11 days suggest that they have the potential to record magmatic processes. Several caveats such as poorly constrained initial conditions and analytical challenges limit the practicality of enrichment peak geospeedometry. Only broad zoning of non-coupled Li should be modeled for timescales of magmatic processes.


INTRODUCTION
Lithium is an increasingly popular element for investigating diverse processes in volcanic systems, including volatile partitioning and degassing (e.g., Vlastélic et al., 2011;Edmonds, 2015), magmatic differentiation (Weyer and Seitz, 2012), and diffusive re-equilibration (e.g., Beck et al., 2006;Parkinson et al., 2007;Gallagher and Elliott, 2009;Ellis et al., 2018). As a fast-diffusing cation, Li re-equilibrates faster than all major and minor elements in olivine and is a potentially powerful geospeedometer for characterizing the sub-solidus histories of magma storage , rapidly occurring volcanic processes (Charlier et al., 2012;Lynn et al., 2018), and the post-eruptive cooling of volcanic deposits (Ellis et al., 2018). However, trace element crystal chemical effects such as charge balancing requirements and element coupling may obfuscate Li geospeedometry. Coupling with other elements leads to multi-component diffusion, where the mobility of an element is also dependent on the concentration and mobility of other elements (e.g., Zhang, 2010).
Lithium diffusion in zircon exemplifies these complexities despite its usage to infer the timescales of magmatic processes in silicic magmas (e.g., Ushikubo et al., 2008;Cooper et al., 2017;Tang et al., 2017;Wilson et al., 2017). Retention of Li in zircon despite crystal ages spanning 10s of thousands of years  sparked debate over apparently "slow" diffusion when Li is coupled with rare earth elements (REE) vs. "fast" when it is not coupled (Ushikubo et al., 2008;Cooper et al., 2017;Tang et al., 2017;Wilson et al., 2017). Trace element zoning in zircon is also strongly dependent on the 3D crystal structure during growth (Trail et al., 2016), making clear distinctions between non-coupled vs. coupled Li diffusion challenging. As a result, timescales and magmatic process interpretations drawn from Li zoning in zircon are currently vigorously debated Rubin et al., 2017;Wilson et al., 2017;Sliwinski et al., 2018).
Trace element zoning in olivine, the dominant phenocryst phase in most basaltic systems, is also complicated by charge coupling effects. Numerous studies have proposed that Li in olivine is subject to charge balancing requirements of other trace elements that are incorporated during rapid crystal growth (e.g., P, Al, Cr; Woodland et al., 2004;Milman-Barris et al., 2008;Mallmann et al., 2009;Grant and Wood, 2010;Spandler and O'Neill, 2010;Tomascak et al., 2016). Monovalent cations (e.g., Li + , Na + , H + ) are thought to occupy octahedrally coordinated sites to charge balance P 5+ substitution in tetrahedrally coordinated sites (Woodland et al., 2004;Milman-Barris et al., 2008;Mallmann et al., 2009) based on the following reaction: IV Si 4+ + VI Mg 2+ = IV P 5+ + VI Li + , Na + , H + where IV denotes coordination in a tetrahedral site and VI denotes that for the octahedral site. Phosphorus, a normally highly incompatible element, is enriched in a skeletal olivine framework during early rapid growth (Milman-Barris et al., 2008;Welsch et al., 2013Welsch et al., , 2014Shea et al., 2019). The above mechanism (Eq. 1) would allow rapidly growing olivine to incorporate monovalent cations like Li + . Phosphorus is also one of the slowest diffusing elements in olivine (Watson et al., 2015), meaning that the mobility of faster diffusing Li (Dohmen et al., 2010) may be subsequently hindered by the charge coupling conditions in Eq. 1. Alternatively, Li may exhibit "inter-site reaction" diffusion behavior analogous to H + decorating preexisting defects associated with trace elements like Ti (Jollands et al., 2016). In this case, Li would diffuse through P-rich regions by hopping between site vacancies or interstitial sites and vacancies associated with trace element charge balancing. Phosphorus enrichments are probably also balanced by a combination of other substitution reactions, including vacancies and/or 3 + ions in octahedral sites (Agrell et al., 1998;Boesenberg et al., 2004;Milman-Barris et al., 2008;Mallmann et al., 2009;Boesenberg and Hewins, 2010): 2 IV Si 4+ = IV P 5+ + IV Al 3+ (4) where VI M 2+ refers to a divalent cation in an octahedral site (e.g., Fe 2+ , Mg 2+ ), VI R 3+ refers to a trivalent cation in an octahedral site (e.g., Al 3+ , Cr 3+ , Fe 3+ ), and VI [] represents an octahedral site vacancy. Aluminum also frequently shows fine-scale zoning inferred to preserve rapid crystal growth and may partially charge balance P 5+ via tetrahedral or octahedral substitutions outlined in Eqs 3 and/or 4 (Evans et al., 2008;Milman-Barris et al., 2008;Welsch et al., 2013Welsch et al., , 2014Zhukova et al., 2017;Shea et al., 2019).
Here, we document the distribution of Li, P, and Al in magmatic olivine to distinguish the signatures of growth-and diffusion-induced zoning for these trace elements in natural basalts. Pumice from explosive Keanakāko'i Tephra (KT, 1500-1820 C.E.; Swanson et al., 2012) eruptions at Kīlauea Volcano are ideal for this study because of the following reasons: (a) the tephra contains euhedral olivine that was easily extracted and oriented perpendicular to a crystallographic axis for sectioning; (b) samples were quenched very rapidly, avoiding complications of post-eruptive Li diffusion (e.g., Ellis et al., 2018); and (c) the magmatic histories of many of the tephra units are well characterized based on previous studies of glass and olivine compositions (Lynn et al., 2017Garcia et al., 2018). We combine backscattered images and X-ray element maps of P and Al by EPMA; laser ablation ICPMS measurements of Li, P, and Al; and nanoscale secondary ion mass spectrometry (nanoSIMS) measurements of Li and P in linear core-to-rim transects and ion element maps and to document the distribution and correlation of Li and P. These data are used to identify coupled vs. noncoupled Li zoning and provide recommendations on selecting Li profiles for geospeedometry.

Sample Preparation and Previous Work
Samples used in this study were previously examined for major and minor elements in olivine (see Lynn et al., 2017Lynn et al., , 2018 and major, minor, and trace element geochemistry in glasses . The euhedral olivine phenocrysts (>0.5 mm) examined in this study had core forsterite [Fo mol%; (Mg/(Mg + Fe) × 100)] contents of 83.9-88.5, with one crystal at 79.2 (Supplementary Material). Each crystal was oriented perpendicular to a crystallographic axis, confirmed via electron backscatter diffraction (EBSD) analyses (see below) and/or recognizable euhedral crystal faces (Lynn et al., 2017). Well-oriented crystals allow analytical profiles to be made parallel to principal crystallographic axes and minimize oblique sectioning through complex 3D distributions of elements like P (Welsch et al., 2014).

Electron Probe Micro Analysis (EPMA)
X-ray maps of Al and P in olivine were obtained using the JEOL JXA-8500F Hyperprobe at the University of Hawai'i. Maps were collected using a 15-kV accelerating voltage and a 300-to 500-nA beam current, with dwell times between 100 and 200 ms/pixel.

Electron Backscatter Diffraction
Crystal orientations were confirmed using the HKL Nordlys EBSD detector on the JEOL JSM-5900LV scanning electron microscope at the University of Hawai'i (UH). Final polishing with a 50-to 70-nm colloidal silica suspension for 4 h on a vibratory polisher improved the quality of electron backscatter patterns (EBSPs). Measurements were made using a 70 • sample tilt, 25-kV accelerating voltage, and a working distance of 16-17 mm. Grids of 10 × 10 individual analysis spaced at least 10 µm apart achieved average mean angular deviation values of <1 • . The EBSPs were processed using the HKL Technology Channel 5 software package, and axis locations were output into lower hemisphere stereographic projections.

Laser Ablation Inductively Coupled Plasma Mass Spectrometry
Trace element analyses for 10 transects in 7 olivine crystals were acquired by laser ablation inductively coupled plasma mass spectrometry (LA-ICPMS) using the same methodology of Lynn et al. (2018). Only seven elements were analyzed in our routine ( 7 Li, 23 Na, 27 Al, 29 Si, 31 P, 52 Cr, and 60 Ni), which allowed significantly longer dwell times on Li (65 ms) than typical for multi-element (e.g., 20+) routines (<10 ms per element). The Ni abundances were measured here to compare with previous EPMA analyses (Lynn et al., 2017) and are not discussed further. Both Na and Cr showed inconsistent relationships with Li and P zoning and are not discussed here. Measurements were made with a Photon Machines Analyte G2 excimer laser connected to an iCAP Q ICPMS in an Ar-He atmosphere at the Earth Observatory of Singapore (Nanyang Technological University). Machine settings were optimized to reduce oxide production and interferences by achieving a ThO + signal intensity of <0.6% of the Th + signal. The laser was run at 70% output, 5.9 J cm −2 fluence, and frequency of 8 Hz with constant voltage. Twenty seconds of background signal (laser off) was acquired followed by 40 s sample ablation.
Measurements utilized 10 × 50-µm spots oriented with the 50µm length parallel to the crystal-melt boundary (Supplementary Figure S1). An 11-µm spacing was used to avoid overlap and to create a 1-µm separation between individual spot analyses. Data reduction in Iolite© (Paton et al., 2011) removed the first and last 3 s of counts for each analysis before converting counts to concentrations. Analyses with fluctuations or spikes in the counts that might represent ablation of Cr-spinel, glass, or fluid inclusions were removed from the dataset before processing. The internal reference isotope 29 Si was used to convert raw counts to concentrations based on the EPMA measurements of SiO 2 in olivine (data from Lynn et al., 2017). The NIST 612 silicate glass was used to calibrate relative sensitivities for Li (40.2 ppm), P (46.6 ppm), and Al (10746 ppm). NIST 610 and USGS BCR-2G were run as unknowns (see Supplementary Material). The three reference materials were analyzed in between each olivine transect to monitor instrumental drift. Analytical precision (2σ) values for repeated measurements of NIST 612 glass are 0.83 ppm for Li (compared with 1.3 ppm given error on preferred GeoReM value), 206 ppm for Al (GeoReM error 212 ppm), and 6.3 ppm for P (GeoReM error 6.9 ppm; Supplementary Material). Measured values were compared with the GeoReM database preferred values (Jochum et al., 2011(Jochum et al., , 2016 to assess data quality. Average 2σ errors for Li measurements of KT olivine crystals were calculated by Iolite© (Paton et al., 2011) and are typically ±0.08 ppm for Li abundances of 1-3, ±6 ppm for Al abundances of 100s of ppm, and ±10 ppm for P abundances of 10s to 100s of ppm.

Nanoscale Secondary Ion Mass Spectrometry
High-resolution nanoSIMS analyses of Li and P zoning were collected for two crystals based on coupled increases of Li and P concentrations observed in LA-ICPMS profiles. The isotopes of 31 P and 7 Li were measured using different beam sources on the Cameca Ametek nanoSIMS 50L at Arizona State University. Samples were sputtered with a 133 Cs + primary beam of 100 pA to measure P and a 16 O − primary beam of 100 pA to measure Li. An electron gun was used to compensate sample charge. Both trace elements were normalized to 30 Si, measured in each routine, in order to monitor possible instrumental drift. Point analyses were collected to identify zones of interest before high-resolution line scans were made. Areas 50 × 50 µm (256 × 256 px) were pre-sputtered to acquire ion images that were used to precisely orient line scans perpendicular to bands of P and Li enrichment. Line scans were collected using a 0.5-µm beam (3 × 3 px) for 500 cycles with a dwell time of 10-20 ms/px and are reported as accumulated counts. Line scans were run past the crystal-melt boundary into regions of adhering glass to isolate mixed signals and remove possible edge effects. The San Carlos olivine standard (NMNH 111312-42; a split specifically characterized for its trace element composition; E. Hauri, personal communication, 2017) was measured with a 5 × 5-µm raster (128 × 128 px) with a dwell time of 1 ms/px for 10 blocks and 10 cycles/block. Phosphorus (23.4 ± 1.2 ppm) and Li concentrations (1.63 ± 0.08 ppm; E. Hauri, unpublished) were used to convert count ratios of 31 P/ 30 Si and 7 Li/ 30 Si to concentrations. There is excellent correlation between the 31 P/ 30 Si ratio and 31 P counts (R 2 = 0.9762) and the 7 Li/ 30 Si ratio and 7 Li counts (R 2 = 0.9998), demonstrating that the ratios depend primarily on their trace element concentrations and were not strongly affected by 30 Si variations or instrument drift (Supplementary Figure S2).

RESULTS
Lithium and P zoning in the Kīlauea olivine crystals follow two systematic and distinct behaviors: non-coupled and coupled.
Non-coupled behavior usually shows a smooth gradient in Li with 1-2 ppm variation over 10s to 100s of µm between the crystal core and rim ( Figure 1A), which is amenable to diffusion modeling. Diffuse, non-coupled Li zoning always occurs in transects where P concentrations are low (a few 10s of ppm) and homogeneous within LA-ICPMS error (typically ∼10 ppm; Figure 1A). We chose to use the term "non-coupled" because it indicates simply that the P and Li concentrations are not correlated. The term "uncoupled" (used by Tang et al., 2017) or "decoupled" implies that the trace elements may have been coupled originally. Non-coupled profiles occur in olivine crystals with both normal (Li core > Li rim ; Figure 1) and reverse zoning (Li core < Li rim ). These profiles from Kīlauea olivine were previously interpreted to represent diffusive re-equilibration initiated by late-stage magmatic intrusions and were modeled to extract the timing and duration of magma mixing and transport .
In "coupled" Li and P zoning profiles, abrupt changes in Li concentration are positively correlated with variations in P ( Figure 1B). Zoning of Li typically increases on the order of 15-80% (0.3-1.6 ppm), compared with P variations that can be up to 2,800% (10s to 100s of ppm) across an entire profile. Six of the analytical profiles (4 of 7 crystals) have coupled Li + P zoning. Coupled Li does not systematically correlate with Al.
Both non-coupled and coupled behaviors of Li + P are observed in some olivine crystals (Figure 2). The core-to-rim transects in these crystals have broad regions of diffuse Li zoning (occurring at P concentrations of only a few 10s of ppm) and FIGURE 1 | Examples of non-coupled and coupled Li behaviors in Kīlauea olivine. (A) Non-coupled behavior of Li. Diffuse Li is strongly zoned, whereas P concentrations are low and essentially non-zoned outside the error bars. The gray line through the Li data is the 1D diffusion model best fit of 18 days, and the dashed lines represent model uncertainties of -4 and +5 days  that reflect the 2σ analytical error and spacing of analyses (1 spot per 11 µm) Additional timescale uncertainties related to temperature calculations (±10 • C; Helz and Thornber, 1987) are on average 10% of the timescales calculated for a given profile. The BSE image (lower left) and EBSD lower hemisphere stereonet projection show that the traverse was done parallel to the a-axis. (B) Coupled behavior of Li. Regions of high P marked by gray bars also correspond to regions of high Li, produced by elemental charge balancing requirements imposed during crystal growth. The orientation of this crystal was estimated by comparing the 2D section with an ideal olivine section taken perpendicular to the c-axis (e.g., Shea et al., 2015a). The traverse was collected roughly parallel to the b-axis.
Frontiers in Earth Science | www.frontiersin.org FIGURE 2 | Two transects that show both non-coupled and coupled behaviors of Li and P. The two profiles (marked by red arrows on BSE image) from the same olivine crystal parallel to the b-and c-axes have diffuse zoning of lithium that can be fit with 1D diffusion models. The gray line through the Li data is the model best fit, and the dashed lines represent model uncertainties that reflect the 2σ analytical error and spacing of analyses (1 spot per 11 µm; Lynn et al., 2018). Additional timescale uncertainties related to temperature calculations (±10 • C; Helz and Thornber, 1987) are on average 10% of the timescales calculated for a given profile. These profiles also contain enrichment peaks (marked by gray bars) that have correlated Li and P enrichments. In this example, Al is also positively correlated with Li and P. The EPMA X-ray maps for Al and P show correlations in enrichment peak regions. The red line on the P X-ray map indicates the location of the nanoSIMS profile shown in Figure 3, which measured the enrichment peak outlined in the dashed black box. The nanoSIMS ion image shows that the enrichment peak (∼15 µm wide in blue on nanoSIMS image) is straddled by two laser ablation pits, demonstrating how small enrichment features can be distorted or homogenized by coarser beam size analytical routines. LA-ICPMS pits (dotted white boxes) appear to have high counts because the ablated crystal in those regions is poorly polished (e.g., not a true reflection of counts or concentration). Analytical error is smaller than the symbol size.
10-to 30-µm-wide regions of elevated Li + P. In such complex profiles, coupled Li + P usually occurs as abrupt increases in concentration relative to adjacent portions of the crystal. These locally coupled regions are herein referred to as "enrichment peaks." LA-ICPMS profiles reveal peaks with a 10-20% (0.1-0.6 ppm) increase of Li and a 300-900% (10s to 100s of ppm) increase of P compared with adjacent analyses (Figure 2). Where present, Li enrichment peaks are always correlated with P enrichments. Aluminum, another trace element in olivine that is commonly but not always positively correlated with P, is not systematically correlated with Li enrichment peaks.
Measurements by nanoSIMS show that enrichment peaks have much higher concentrations (e.g., 550 vs. 340 ppm; b-axis profile in Figure 2 vs. Figure 3) than measured by LA-ICPMS because of the much smaller analytical spot size (e.g., <1 µm vs. 10 µm, respectively). The shape of this peak is better defined in the nanoSIMS profile, which shows a doublet in P zoning (Figure 3) that is not resolved by LA-ICPMS (Figure 2). Another enrichment peak that is 15 µm wide in a Unit K1 olivine has 700 ppm P measured by nanoSIMS, but the peak is absent in the LA-ICPMS profile from the same crystal (Figure 4).
The coupled vs. non-coupled behaviors can be distinguished by distinct trends in P vs. Li expressed as cations per four oxygens (atoms per formula unit; Figure 5). Non-coupled behavior shows a near-vertical trend with variable Li at relatively constant P. Coupled Li and P covary positively with a shallow slope and relatively high R 2 values ( Figure 5). Analytical profiles that exhibit both behaviors can be divided into two groups of data with either a near-vertical trend or a shallow positive trend ( Figure 5C).

DISCUSSION
Coupled and non-coupled behaviors of Li with P zoning in olivine clearly distinguish signatures of crystal growth and diffusive reequilibration. These features provide a framework for applying Li geospeedometry in olivine. The abundances, zoning styles, and spatial variations in Li and P correlations are used in the following sections to (a) examine the value of high-resolution analytical methods for discerning and interpreting small-scale zoning features to understand trace element behavior in minerals; (b) provide insights into the mechanisms controlling non-coupled and coupled Li zoning; (c) illustrate the complex, 3D distribution of these elements; (d) test diffusion models on coupled Li and P enrichment peaks; and (e) provide guidelines for applying and interpreting Li geospeedometry in olivine.

High-Resolution Is Needed to Characterize Enrichment Peak Features
The spatial resolution of trace element analyses affects the magnitude of apparent enrichment relative to non-zoned portions of crystals (e.g., Bouvet de Maisonneuve et al., 2016;Shea et al., 2019). Trace element enrichment peaks in olivine that reflect rapid crystal growth are typically less than a few 10s of µm wide (e.g., Welsch et al., 2013Welsch et al., , 2014Shea et al., 2015bShea et al., , 2019Manzini et al., 2017). Analytical methods with large beam sizes (e.g., this LA-ICPMS method with 10-µm-wide spots) or continuous scanning (e.g., raster method; Bouvet de Maisonneuve et al., 2016) smear fine-scale trace element zoning and yield reduced concentrations in these peaks. Finer beam size methods such as nanoSIMS (<1 µm beam size) allow details of the enrichment peaks to be revealed (Figures 2, 3; see also Manzini et al., 2017). Two olivine crystals were analyzed with LA-ICPMS and nanoSIMS methods to assess the loss of information that occurs when using coarse analytical methods such as laser ablation. The LA-ICPMS profile shows a 20-µm-wide Li + P enrichment peak near the crystal rim with P concentrations of 340 ppm (Figure 2), enriched by about a factor of ∼7 compared with the nearby non-zoned regions characterized by ∼50 ppm P. The nanoSIMS analyses of the same peak revealed a P doublet ∼10 µm wide with up to 550 ppm (Figure 3). The 60% higher concentration of P and presence of a doublet in the nanoSIMS profile illustrate the superiority of this method for resolving fine-scale details about trace element zoning in olivine as well as other minerals (e.g., Manzini et al., 2017;Tang et al., 2017;Seitz et al., 2018). Thus, coarse-resolution methods such as laser ablation are not suitable for characterizing and subsequently modeling the enrichment peak features.

Diffusion Generates Broad Non-coupled Lithium Zoning
The majority of Li profiles are dominantly non-coupled and amenable to diffusion modeling. In these profiles, P concentrations are low, a few 10s of ppm, and there is no correlation of P with Li content (Figure 5). Thus, we propose here that Li is not charge balancing P or other elements at such The BSE image (upper left) and EBSD lower hemisphere stereonet projection show that this crystal is elongated along the a-axis and was sectioned perpendicular to the c-axis. (B) Overall diffuse zoning of Li and low P concentrations, as measured by LA-ICPMS with modeled timescale and uncertainties associated with analytical error and spatial resolution (described above). (C,D) NanoSIMS ion images show that proximal to the olivine-melt boundary, there is a <10-µm enrichment peak with coupled Li and P. (E) NanoSIMS measurements across the enrichment peak, showing coupled P and Li in the peak region and steadily increasing, non-coupled Li surrounding it. Analytical error on all nanoSIMS data is smaller than the symbol size.
low concentrations (typically <50 ppm) and is able to diffuse readily. Charge balancing for P substitution in the tetrahedral site is balanced by other ionic species (Al, Cr, Fe 3+ ) and/or vacancies (Eqs 2-4). Aluminum charge balances most or all of the P budget in low-P regions of both natural and experimental olivine . Regions of both low P and Al (e.g., <100 molar ppm) have an approximately 1:1 correlation (region "3" on Figure 10b in Shea et al., 2019). Thus, P is balanced by Al in low-P regions, and Li is unaffected by charge balance requirements.
In the absence of charge balancing requirements, we propose that Li atoms diffuse within the olivine lattice via the slower vacancy diffusion mechanism of 7 Li inferred for low-Li basaltic systems (Dohmen et al., 2010;Lynn et al., 2018). This interpretation is supported by our previous documentation of both normally and reversely zoned Li profiles in diverse olivine populations within these samples . This indicates a common origin for the non-coupled Li profiles via reequilibration after magma mixing, yielding timescales that range from hours to days that were interpreted to reflect late-stage magma mixing and ascent.

Charge Coupling of Li and P in Enrichment Peaks Record Crystal Growth
Lithium enrichment peaks are always positively correlated with P enrichments that are 10s to 100s of ppm higher than in regions of low P (e. g., Figures 2, 4). This coupled behavior is consistent with the charge balancing substitution reaction in Eq. 1 and therefore reflects growth, not diffusion. The high-resolution nanoSIMS measurements show systematic and positive correlation between P and Li within an enrichment peak (Figure 5C), indicating that Li in these regions is likely utilized for charge balancing. The similar widths of the Li and P peaks support this interpretation ( Figure 4E). Lithium is also two orders of magnitude less abundant than P measured in the enrichment peaks (1-3 ppm vs. 100s of ppm; Figure 5), indicating that charge coupling will dominate the behavior of Li in these regions. This is consistent with interpretations of ppb level Li contents charge coupling with 10 2 -10 4 ppm concentrations of Y + REE in zircon Tang et al., 2017;Wilson et al., 2017). Thus, we interpret that Li enrichment peaks are generated during crystal growth as Li + ions in octahedrally coordinated Mg sites in olivine charge balance P 5+ ions occupying tetrahedral sites (Woodland et al., 2004;Mallmann et al., 2009).
Lithium enrichments are inconsistently correlated with Al, another trace element that is commonly but not always associated with P zoning (Evans et al., 2008;Milman-Barris et al., 2008;Welsch et al., 2013Welsch et al., , 2014Zhukova et al., 2017;Shea et al., 2019). Uncommon correlations of Li and Al are likely a consequence of multiple charge balancing relationships that share P as a common element (e.g., Eqs 1, 3, and 4). Aluminum charge balances most of the P budget, the behavior of Al during crystal growth is different from P, and Al is frequently a marker for sector zoning in olivine, contrary to P zoning (Pack and Palme, 2003;Milman-Barris et al., 2008;Shea et al., 2019). Thus, the non-systematic correlation of Li with Al in enrichment peaks probably results from several charge balancing equations accommodating most of the budget of P-rich zones.

A Model of Coexisting Growth and Diffusion Signatures
Interpreting non-coupled vs. coupled Li signatures in olivineand subsequently extracting timescales of magmatic processes from Li zoning-requires careful consideration of zoning variations in individual transects and distributions across 2D olivine sections. Current models of P zoning in olivine (e.g., Welsch et al., 2014) demonstrate how coupled and non-coupled (i.e., growth and diffusion) Li zoning patterns may occur within individual transects and single olivine crystals (Figure 6).
The initial dendritic or hopper morphology of rapidly growing olivine (Donaldson, 1976) is recorded by P-rich (and by extension, Li-rich) zones ( Figure 6A). Phosphorus zoning is concentrated during rapid crystal growth in both primary and secondary branches that are usually <10 µm wide; they represent a minor fraction of the olivine volume ( Figure 6B, following Welsch et al., 2014). Most of the olivine volume has low P, as FIGURE 6 | Illustration of Li and P zoning resulting from growth and diffusion. (A) Euhedral olivine crystal with volumetrically small P-rich skeletal structure (e.g., Welsch et al., 2014). Dashed line marks the plane of the section perpendicular to the b-axis. (B) Zoning of Li + P rich skeleton (light blue) from rapid growth surrounded by Li-and P-poor compositions from later slow growth during crystal maturation (dark blue). Li-and P-rich zones in secondary branches may be discontinuous or continuous. White lines mark the locations of compositional profiles presented below the section. (C) Zoning generated by subsequent diffusion of Li (pink rim). Transect 1 has a diffusive rim with no change to the Li enrichment peak. Transect 2 shows Li diffusion in the absence of growth-related enrichment peak zoning. P is unchanged throughout the diffusion process due to its slow diffusivity (Watson et al., 2015) and short time elapsed (hours to days) in this scenario. maturation of an initially skeletal crystal leads to in-filling under slower growth conditions. Thus, most of the area in 2D sections will have low P concentrations (a few 10s of ppm; Figure 6B). Because Li is interpreted to not charge balance low-P regions, non-coupled Li throughout most of the olivine volume permits diffusive re-equilibration as a dominant signature for zoning ( Figure 6C, a-axis profile). Due to the complex distribution of P zoning in 3D, analytical traverses within a single 2D section may sample both fine-scale P-rich zones and low-P regions (Figure 6C, c-axis profile). Thus, individual profiles that exhibit both non-coupled and coupled Li behaviors record both crystal growth and subsequent diffusion processes.
Multiple transects collected within a single olivine crystal in this study record different Li and P zoning patterns as a result of this complex 3D distribution of trace elements. Within the olivine crystal in Figure 4, LA-ICPMS and nanoSIMS profiles located 100s of µm apart from each other record different trace element zoning patterns. Lithium exhibits dominantly non-coupled behavior in the LA-ICPMS profile where maximum P concentrations are up to ∼80 ppm. In contrast, the nanoSIMS profile reveals a clear enrichment peak of 700 ppm P and 3 ppm Li about 8 µm away from the crystal rim ( Figure 4E). Differences in zoning between these profiles may be due to the loss of information from the 10µm-wide LA-ICPMS analyses traverse (as discussed above). However, 2D sectioning through complex 3D distributions of trace elements likely also plays a role. The nanoSIMS ion image of P zoning in Figure 4E shows clearly decreasing counts along the length of the secondary branch toward the direction of the LA-ICPMS profile ( Figure 4A). Thus, the Li-and P-rich zone sampled by the nanoSIMS traverse in Figure 4 is likely not continuous along the crystal face and was not sampled by the LA-ICPMS traverse. This example underscores the complexities of interpreting 1D transects that sample complex 3D distributions of trace elements in olivine. Collecting multiple traverses from well-oriented crystals is essential for using trace element zoning for robust interpretations of geologic processes.

Diffusion Models of Enrichment Peaks
Transects exhibiting both coupled and non-coupled Li behaviors (Figure 2) contradict the notion that rapid diffusion of Li in olivine should homogenize zoning over hours to days within basaltic systems. Here we apply diffusion models to the highly resolved nanoSIMS enrichment peaks of Li and P, as has been done for Li in zircon (e.g., Rubin et al., 2017), to assess the potential impacts of coupled diffusion on Li timescales. The results are used to determine whether these fine-scale features can be used to extract meaningful information about magmatic processes.
Coupled Li in enrichment peaks is likely controlled by the slower diffusion of P (Watson et al., 2015) because charge balanced elements within the olivine lattice are susceptible to sharing coupled fluxes during diffusive re-equilibration (Lasaga, 1979;Ganguly, 2002). Several diffusion models for the nanoSIMS Li profiles were run assuming either non-coupled or coupled behavior to test whether coupled fluxes influence Li diffusion in enrichment peaks. Non-coupled Li diffusion proceeds at the expected "slow" rate for Li in basaltic systems following the vacancy mechanism of Dohmen et al. (2010). Coupled Li diffusion is assumed to be slowed to the rate of phosphorus diffusion (D P ; Watson et al., 2015) due to charge coupling. These assumptions yield radically different diffusion timescales because the diffusion coefficient for phosphorus at 1,200 • C is 3.5 orders of magnitude slower than Li (D P = 10 −18.1 m 2 /s vs. D 7Li = 10 −14.6 m 2 ; Watson et al., 2015;Dohmen et al., 2010).
Finite difference diffusion models were applied to the coupled nanoSIMS peaks in Figure 4 to evaluate whether one or both of these diffusion coefficients are viable for applying Li geospeedometry in enrichment peaks. We used the nonconcentration dependent one-dimensional form of Fick's second law (Crank, 1975): where C is concentration (ppm), t is time (s), D is the diffusion coefficient (m 2 /s), and x is distance (m). The diffusion coefficient of Li (D 7Li ) parallel to the c-axis can be calculated using the Arrhenius relation (Dohmen et al., 2010): where T is temperature. This expression is appropriate for diffusion in olivine containing 1-10 ppm Li and temperatures 800-1,200 • C (Dohmen et al., 2010). The diffusion coefficient of P in olivine can be calculated using the Arrhenius relation (Watson et al., 2015): where R is the ideal gas constant (8.314 J/mol K). The Li enrichment peak was modeled using both D 7Li and D P to assess potential timescales of non-coupled and coupled diffusion, respectively. The model temperature was assumed to be constant at 1,237 • C (taken from Li diffusion models in Lynn et al., 2018) and was calculated using the Kīlauea glass MgO thermometer (uncertainty of ±10 • C; Helz and Thornber, 1987) and the inferred melt Mg# in equilibrium with olivine rim Fo contents (Kd ol/melt Fe−Mg =0.343 ± 0.018; Matzen et al., 2011). This value is akin to the magma temperature for mixing and diffusion recorded in the near-rim olivine composition (Lynn et al., 2017. This approach is suitable for our modeling here because the nanoSIMS enrichment peak is located very near to the crystalmelt boundary (Figure 4).
Initial conditions (C 0 ; the composition of the enrichment peak prior to diffusion) have a high degree of uncertainty, and several were estimated due to the absence of a compositional plateau in the enrichment peak. We assessed this effect on retrieved timescales by testing four assumptions: (1) the initial peak had a flat plateau near the maximum measured concentration, (2) the initial peak had a "dipping plateau" akin to sectioning the enrichment peak off-center from its maximum value, (3) the initial peak plateau was higher than maximum measured value, and (4) the initial peak was tapered with a maximum value higher than that measured in the profile ( Table 1). An additional model was run to fit the non-coupled region of the profile assuming that Li ions diffused through the enrichment peak without significantly modifying it. This follows the hypothesis of inter-site reaction diffusion, where Li ions may diffuse quickly through P-rich regions by hopping between site vacancies and/or interstitial sites and vacancies associated with P. This is analogous to previous studies of H + diffusion in olivine, where H + was decorating pre-existing defects associated with trace elements like Ti (Jollands et al., 2016). Boundary conditions were set to P and Li concentrations in the far field around the enrichment peaks, and Li values were constrained by the core composition measured by LA-ICPMS in Figure 4 and Table 1. 1 | Timescales of diffusive re-equilibration for diffusion models using the Unit K1 olivine 3 coupled enrichment peak in Figure 4.

Interpreting Enrichment Peak Timescales
Modeling the Li enrichment peak in Figure 4 using D Li (Dohmen et al., 2010) yields a diffusion timescale of 2.8-17 min depending on the initial condition assumptions (Table 1 and Figure 7). These timescales are far too short to reflect a magma mixing event prior to eruption, and they do not agree with timescales from non-coupled Li profiles (hours to days; Lynn et al., 2018). It is also unlikely that they represent the time elapsed between crystal growth and eruption. The short Li timescales are also inconsistent with diffusion timescales recorded by the P peak (using D P ; Watson et al., 2015), which are much longer and range from 5.5 to 13 days ( Table 1). Inconsistent Li and P timescales modeled here support the interpretation that Li + P are likely coupled in enrichment peak regions. Thus, Li enrichment peaks should not be modeled using D Li . If instead the Li peak is fit using the diffusion coefficient for P (Watson et al., 2015), the resulting Li timescales are broadly congruent with P and range from 5.0 to 28 days ( Table 1). Regardless of the assumed initial conditions, congruent timescales derived from using D P to model both profiles suggest that Li in enrichment peaks is subject to coupled diffusive fluxes with P. These durations are also generally consistent with the latestage magma mixing timescales of hours to days obtained from non-coupled Li zoning in other Kīlauea olivine crystals (generally hours to days, but up to a few weeks; Lynn et al., 2018). Thus, diffusion chronometry of enrichment peak features in olivine might need to account for coupled diffusive fluxes by using D P to model coupled Li behavior.
An additional model was run to fit the non-coupled portion of the nanoSIMS profile assuming that Li ions diffuse through the enrichment peak without significantly modifying it. Following this inter-site reaction hypothesis, the enrichment peak does not hinder the movement of Li ions, and the timescale retrieved reflects the non-coupled diffusion of Li from the olivine core to its rim. The retrieved timescale is 9.6 h, but we note that the 25µm nanoSIMS profile lacks a clear core plateau and only samples a small part of the olivine rim. The 9.6-h duration is shorter than the timescale of several days modeled for the full crystal's LA-ICP-MS profile in Figure 5, which is likely a symptom of incompletely measuring the full diffusion profile from core to rim.

Caveats and Challenges With Modeling Enrichment Peaks
Several caveats and challenges limit the application of enrichment peak chronometry despite the congruency between Li and P timescales derived from our nanoSIMS profiles. Interpreting timescales from enrichment peaks is made difficult by their 3D distribution within the olivine volume during initial rapid growth (e.g., Welsch et al., 2014). The enrichment peak from Figure 4 is located within 15 µm of the crystal-melt boundary and could plausibly be related to a late-stage process occurring shortly before eruption. However, there is increasing consensus that mineral chemistry from core to rim (following the canonical "tree-ring" approach) does not provide chronological records of the magmatic histories (e.g., Welsch et al., 2014;Shea et al., 2015bShea et al., , 2019. Phosphorus-rich zones located near crystal rims can be products of early stages of crystal growth, invalidating the geologic interpretation of the timescales retrieved here. Model temperatures are also uncertain, as the thermal history associated with the early crystal growth is not well recorded and the entire crystal history is likely not isothermal. Temperatures calculated for Kīlauea magmas using the glass MgO thermometer (Helz and Thornber, 1987) and the inferred melt MgO in equilibrium with the olivine forsterite content (using the Fe-Mg partition coefficient K D = 0.343 ± 0.008 of Matzen et al., 2011) likely record a more recent temperature history for the olivine crystals than the early P-rich zones experienced. Thus, a more detailed examination of P zoning in high resolution, the timescales derived from enrichment peaks, and their distribution across 2D sections is needed before using these features to characterize magmatic histories.
The enrichment peak diffusion timescales are also highly uncertain due to the absence of initial plateaus in our nanoSIMS data. Our retrieved timescales are largely insensitive to the shape of our initial conditions (flat vs. dipping plateau) but are strongly affected by the concentration chosen for the plateau (Figure 7). At present, few studies have measured P zoning profiles in high resolution (Bouvet de Maisonneuve et al., 2016;Manzini et al., 2017), and those profiles generally lack clear compositional plateaus. Fewer P profiles have been modeled for diffusive reequilibration (e.g., Manzini et al., 2017), and uncertainty related to unknown initial conditions was not addressed. Because the P-rich zones introduce an extremely steep compositional contrast with nearby low-P regions (100s vs. 10s of ppm), P may never preserve an initial plateau, and constraining an appropriate initial condition would be impractical. The assumption of a step function across branch may also be inappropriate-growth models of fine-scale P-rich zones show they are subject to micrometer-scale P zoning at their edges instead of sharp boundaries with adjacent low-P regions . In this FIGURE 7 | Diffusion models of Li and P enrichment peaks with different assumed initial conditions (C 0 ; dashed lines). Analytical error on all nanoSIMS data is smaller than the symbol size. Model fits are shown in solid gray lines. Lithium peaks were modeled assuming that the sloped region surrounding the enrichment peak was generated by non-coupled Li diffusion and is not fit here. Initial conditions were constrained by the Li content at the core of the olivine crystal (1.5 ppm far field; from LA-ICPMS data in Figure 5) and different assumptions about peak compositions: (A) a flat plateau near the maximum value measured in the profile (C 0 1), (B) a "sloping" plateau from a non-ideal section through the enrichment peak (C 0 2), (C) a higher plateau than measured (actual value unknown; C 0 3), and (D) a tapered initial condition with a maximum value greater than was measured (actual value unknown; C 0 4). Li timescales represent models that utilized the diffusion coefficient for Li (t Li ; Dohmen et al., 2010) and for P (t P ; Watson et al., 2015). (E) Additional model to fit the non-coupled portion of the profile assuming inter-site reactions allows Li ions to diffuse through the enrichment peak without modifying it significantly (after models for H diffusion; Jollands et al., 2019). case, the measured enrichment peaks may dominantly reflect crystal growth and minimal diffusive re-equilibration.
A remaining challenge with enrichment peaks is understanding how Li ions diffuse within peaks and/or the regions surrounding them. The assumptions of coupled diffusion for the models above follow recent studies on zircon in assuming that Li mobility is dependent on P in these regions and diffusion through the peaks is slowed relative to adjacent regions. However, Li may hop between sites and vacancies in the olivine lattice via an inter-defect site reaction behavior, as has been shown for H + (e.g., Jollands et al., 2016Jollands et al., , 2019. If Li diffusion occurs via site reactions, then the assumption of Li and P diffusive coupling (D Li = D P ) is invalid. Inter-defect site reactions would allow Li to diffuse through the peak regions, rather than being slowed within them. Regardless of how Li behaves in/around peak regions, we maintain that enrichment peak features should not be modeled for timescales.
Finally, collecting enough Li and P enrichment peak profiles for accurate timescales modeling is both analytically challenging and potentially expensive. At least 20 profiles are needed to accurately constrain diffusion timescales in a given olivine population, and complex zoning or magmatic histories warrant more (Shea et al., 2015a). Measurements of Li and P by nanoSIMS require different primary beams ( 16 O − and 133 Cs + , respectively), so the minimum number of timescales suggested by Shea et al. (2015a) requires at least 40 profiles measured by nanoSIMS to characterize 20 enrichment peaks. Given the caveats, challenges, and analytical difficulty of measuring and modeling nanoSIMS profiles of Li and P, fine-scale enrichment peaks in olivine are unlikely to be viable chronometers for most studies.

Guidelines for Applying Lithium Geospeedometry in Olivine
The dual behavior of both coupled and non-coupled Li zoning within individual olivine crystals provides a framework on how to apply Li geospeedometry in basaltic systems. Diffuse, noncoupled profiles of Li should be modeled using D Li (vacancy mechanism from Dohmen et al., 2010;Lynn et al., 2018). This approach is applicable to both normally and reversely zoned lithium profiles to provide timescales for magma mixing events shortly (hours to days) before eruption in basaltic systems. Lithium enrichment peaks in olivine should only be modeled if high-resolution analyses resolve many peaks well (e.g., with nanoSIMS), an initial plateau can be identified, and better constraints on the diffusion mechanism controlling Li in peak regions (e.g., coupled, inter-defect site reactions) can be defined. Ideally, these coupled Li timescales should also be interpreted in context with other diffusion chronometers.
The complexities arising from coupled Li behavior may be avoided by modeling only non-coupled Li profiles that are greater than a few 10s of µm wide, minimizing the risk of increasing Li at the rim caused by coupled enrichments of P and Li. Coupled Li and P can also be avoided by preferentially selecting and measuring well-oriented olivine crystals. Crystals with obvious skeletal morphologies (Figure 6) should be avoided because they represent the blueprint of crystal growth with pervasive P-rich zones. In the case of thin sections, crystals should be screened for those closest to perpendicular to the a-or b-axis, and 1D traverses should be taken in the center of well-defined crystal faces to take advantage of incomplete secondary branching (e.g., Figure 4).
Ideally, euhedral olivine crystals should be sectioned through crystal cores, perpendicular to a crystallographic axis. Preferential sectioning of olivine crystals perpendicular to the a-axis and traverses parallel to the b-axis would be most likely to yield noncoupled Li profiles for modeling. Olivine crystals are commonly initially elongated along the a-axis, and branches rich in Li and P are concentrated in the a-c plane during rapid growth . Subsequent maturation and slower olivine growth occur along the b-axis, resulting in low-P regions that allow diffusive re-equilibration of non-coupled Li. Sectioning perpendicular to the b-axis and measuring parallel to the c-or a-axis may reduce the chances of crossing P-rich zones (Welsch et al., 2014;cf. our Figure 1A).

CONCLUSION
Carefully oriented and sectioned olivine crystals provide a unique opportunity to document trace element zoning and the timescales of processes in basaltic magmas. Profiles using Li zoning are especially important for characterizing magma mixing events just prior to eruption (hours to days). Combining profiles of Li and P zoning using EMPA, LA-ICPMS, and nanoSIMS methods in rapidly quenched olivine from explosive Kīlauea eruptions reveals systematic differences between records of crystal growth and subsequent diffusive re-equilibration.
(1) Lithium and phosphorus have two distinct relationships in olivine: (a) non-coupled, wherein Li diffuses over 10-300 µm of low P concentrations (a few 10s of ppm) and is amenable to diffusion modeling, and (b) coupled, in which Li and P covary positively and reflect a framework of rapid crystal growth.
(2) Both non-coupled and coupled behaviors can be preserved in a single crystal because the dendritic nature of early rapid growth results in trace element enrichments that are volumetrically limited within the olivine. The preservation of Li and P enrichment peaks within diffusion profiles is a consequence of Li maintaining local (e.g., few µm scale) growth-induced charge balancing with P. (3) Li is always systematically coupled to P enrichments, although P zoning may also be correlated with other trace elements (e.g., Al). We interpret this to reflect crystal growth processes and the role of P in driving the charge balancing requirements in the olivine lattice. (4) Robust interpretations of trace element zoning documented by 1D analytical transects require knowledge of the 2D and 3D spatial distribution given the complexity of zoning generated during crystal growth.
These observations provide a framework for identifying Li profiles suitable for diffusion modeling. The following practical guidelines are recommended: • Avoid crystals with skeletal or dendritic morphologies.
They probably contain pervasive P-rich zones with coupled Li. • Euhedral crystals should be preferentially sectioned perpendicular to the a-or b-axis. For randomly sectioned olivine (i.e., thin sections), crystals should be screened to avoid sections highly oblique to the a-or b-axis. • Analytical traverses should be taken parallel to crystallographic axes and in the center of well-formed faces to maximize potential of sampling between incomplete secondary branches. Two optimal orientations for measuring non-coupled Li are (1) perpendicular to the a-axis with a traverse measured parallel to the b-axis and (2) perpendicular to the b-axis with traverses measured parallel to either the a-or c-axis. • Collect Li and P data concurrently to discern if profiles are significantly affected by charge coupling. If P data are unavailable, modeling diffuse non-coupled Li that is both normally and reversely zoned within a single eruption (e.g., as a result of magma mixing) is suitable. • Model only wider non-coupled profiles (many 10s of µm) with a well-defined core plateau delineating the initial conditions.
The complexities documented here for olivine highlight the need for careful geochemical studies using well-oriented crystals to resolve the disparate interpretations drawn from trace element zoning in other silicate minerals such as zircon.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
KL conducted LA-ICPMS, EPMA, and EBSD analyses. KL and MG collected nanoSIMS data. All authors assisted in sample collection, contributed to the interpretations. KL led the writing of the manuscript.

FUNDING
This work was supported by the Fred M. Bullard Graduate Fellowship at the University of Hawai'i and the National Science Foundation (NSF) East Asia and Pacific Summer Institutes grant OISE1513668 to KL. Additional support was provided by NSF grants EAR1347915 and EAR1449744 to MG, and EAR1725321 to TS. This is SOEST contribution number 10977.