Millennial storage of near-Moho magma

See allHide authors and affiliations

Science  19 Jul 2019:
Vol. 365, Issue 6450, pp. 260-264
DOI: 10.1126/science.aax4092

Millennial magma reserves

Lava erupts after being stored as magma deeper underground. Silica-rich magmas characteristic of volcanoes like Pinatubo are stored for thousands to hundreds of thousands of years. However, we have little understanding of the time scales for deeper basaltic magma like that supplying Icelandic volcanoes. Mutch et al. used diffusion of chromium and aluminum atoms to show that these magmas are stored for hundreds to thousands of years at the crust-mantle boundary. This discovery elucidates the time scale for understanding how magma is created and stored and how it erupts.

Science, this issue p. 260


The lower crust plays a critical role in the processing of mantle melts and the triggering of volcanic eruptions by supply of magma from greater depth. Our understanding of the deeper parts of magmatic systems is obscured by overprinting of deep signals by shallow processes. We provide a direct estimate of magma residence time in basaltic systems of the deep crust by studying ultramafic nodules from the Borgarhraun eruption in Iceland. Modeling of chromium–aluminum interdiffusion in spinel crystals provides a record of long-term magmatic storage on the order of 1000 years. This places firm constraints on the total crustal residence time of mantle-derived magmas and has important implications for modeling the growth and evolution of transcrustal magmatic systems.

Understanding the temporal evolution of magmatic systems from micrometer to crustal scales has become a major focus in volcanology (1), driven by important advances in the use of crystal chemistry as a magmatic chronometer (2, 3). Compositional zonation in volcanic crystals has been integrated with diffusion models to show that the time scales of days to years preserved in the crystal record can be linked directly to volcanic monitoring data (47). These short time scales are conceptually linked to magma mixing, remobilization of crystal mush, and magma rise from storage regions before eruption. However, the time scales of magma storage and the behavior of magma plumbing systems during periods of apparent quiescence are not well understood. Estimates of magma and mineral residence times in the crust have ranged from 102 to 105 years (8). Parsing these time scales into different crustal levels has proven difficult but is key to understanding transcrustal magmatic systems. Storage time scales in shallow silicic systems can span the aforementioned range in near-solidus conditions of cold storage (8, 9) or in an eruptible high–melt-fraction state (10). By contrast, little is known about storage times in the mafic systems that dominate magmatism at mid-ocean ridges, ocean islands, and in the deeper part of arc volcanoes. Direct storage times in the lower crust are notoriously difficult to estimate by diffusion chronometry as shallow and mid-crustal processing obscure lower-crustal signals through crystal equilibration and resorption. We estimate magma residence time in the lower crust by applying diffusion chronometry to observations of Cr number (Cr#) [Cr# = Cr/(Al + Cr)] zoning in spinels contained in ultramafic nodules from the Borgarhraun eruption in Iceland. This eruption provides constraints on magmatic processes operating close to the crust-mantle boundary, called the Moho.

Borgarhraun is an early postglacial (10.5 to 7 thousand years before the present) primitive basalt flow from the Theistareykir System in the Northern Volcanic Zone of Iceland (11). Geothermobarometry and melt inclusion studies have shown that the magma crystallized near the Moho (depth ~24 km) from a compositionally diverse set of primary mantle melts (1113). Fresh tephra from this eruption contains primitive crystals of olivine (Fo87–92, where Fo is forsterite content in mol. %), clinopyroxene, Cr spinel, and plagioclase (An80–90, where An is anorthite content in mol. %). Mafic and ultramafic cumulate nodules are abundant with textural characteristics typical of storage and maturation in a crystal mush, such as grain rounding and curved crystal edges (14). Wehrlitic nodules show olivine and Cr spinel crystals enclosed in clinopyroxene (Fig. 1), a feature observed in mafic plutonic rocks from layered intrusions (14), oceanic sections (15), and ophiolites (16). The clinopyroxene-hosted olivine and spinel crystals are associated with small melt pockets (Fig. 1) and are surrounded by thin layers of melt (~2 to 10 μm thick). We measured spinel widths of 20 to 600 μm. Spinels with widths below 200 μm have homogeneous Cr# contents that vary between different nodules. Crystals larger than 200 μm are compositionally zoned with rims similar in composition to the small spinels and cores of higher Cr#. The Cr# of these spinel cores scales with spinel size (Fig. 1). The eradication of compositional heterogeneity in the smaller spinels and the preservation of zonation in the larger spinels are likely to be controlled by a diffusive process. Smaller crystals completely equilibrated with their local environment, whereas larger crystals preserved higher Cr# core compositions in a state of disequilibrium. We observed this effect in a large population of spinels (~100) across seven different nodules, which rules out the role of sectioning effects (figs. S1 and S2).

Fig. 1 Disequilibrium in cumulate nodules from beneath the Icelandic Moho.

(A and C) Aluminum maps of Borgarhraun wehrlitic nodules with poikilitic texture and (B and D) accompanying plots that show spinel core Cr# compositions with respect to their measured width. Al maps are colored based on phase: yellow-gray corresponds to Cr spinel (Sp), red corresponds to clinopyroxene (Cpx), black corresponds to olivine (Ol), and gray corresponds to glass (Gl). The black curves are 3D spherical spinel diffusion models (modeled at 1215°C) conducted for different crystal radii and diffusion time scales (years). Error bars correspond to the standard deviation of multiple compositions sampled from each spinel core. (A) and (B) correspond to sample BORG_NOD1_N3 and (C) and (D) to BORG_NOD3_N2.

We expect the crystallization of the spinels and the surrounding clinopyroxene to take place within 1 year of magma emplacement (17). We used thermodynamic modeling to account for the spinel Cr# compositions by equilibration over a cooling interval of ~75°C (figs. S3 and S4), similar to the crystallization interval of the Borgarhraun olivines (1350 to 1215°C) (18). The spinels were initially encased in clinopyroxene at a high temperature and subsequently equilibrated with the surrounding melt pockets and grain-boundary films following cooling to a lower temperature (1215°C) (fig. S5). We infer very little deformation of the spinels as they retained their equidimensional forms, which precludes stress-directed lattice diffusion of Cr and Al to create the concentric zoning patterns that we observed (19).

We cannot resolve the geometry, extent, and interconnectivity of the melt network from two-dimensional (2D) sections of the nodules. However, previous studies have suggested melt pockets and interconnected grain-boundary liquid films play an important role in the final stages of crystallization of large mafic intrusions (14, 20). The melt networks may also have been more extensive before clinopyroxene crystallization, given the important role of the liquid in mush homogenization (21). The composition (e.g., Mg# and Cr#) of the melt pockets adjacent to the spinel crystals was different from that of the final carrier liquid (fig. S6), which suggests that these melt networks were isolated from the open magma body. The variability in spinel Cr# between different nodules highlights heterogeneity in the mush and the development of microenvironments for spinel equilibration.

We used the FEniCS (22) finite element code to model 2D diffusion of Cr–Al exchange in spinels, combined with a nested sampling Bayesian inversion (23) to constrain the timescales of mush storage with robust estimates of uncertainty. We derived diffusivity equations for Cr–Al exchange in Cr spinel that take into account the covariance between the underlying diffusive parameters (24, 25) (figs. S7 and S8 and table S1). We ran the models using Gaussian prior distributions with 1σ uncertainties for temperature (1215 ± 30°C), pressure (0.8 ± 0.14 GPa) (12, 18), and ferric iron contents (Fe3+/Fetotal) of the system (0.14 ± 0.02) (26). We used a multivariate Gaussian prior distribution for the parameters of the Cr–Al interdiffusion coefficient as constrained by its covariance matrix (table S2). We used the equilibration temperature of spinel rim compositions (1215°C) that corresponded to a cooler crystal mush. We assumed the same initial core Cr# of 0.41 for all crystals on the basis of the largest measurable spinel that was most likely to retain its original core composition (Fig. 2). This provided us with a maximum estimate of the storage time scale.

Fig. 2 Timing diffusive equilibration of Borgarhraun spinel crystals.

(A and C) 1D Cr# profiles measured by electron microprobe microanalysis across Borgarhraun spinels. Blue curves show the modeled 1D fit from the 2D finite element diffusion model (insets). Dashed lines are model initial conditions. The error bars are 1σ standard deviations based on analytical uncertainties. (B and D) Temperature-time density plots for posterior distributions from the Bayesian inversion. The quoted time is the median with 1σ uncertainties. (A) and (B) correspond to profile BORG_NOD1_N3_SP1_P1 and (C) and (D) to BORG_NOD3_N2_SP1_P1.

We applied our Bayesian inversion to seven spinel crystals (figs. S9 to S19 and table S3) from three different nodules that converged to magma storage time scales with a median of 1400 years, with 95% of timescales being less than 4100 years (Fig. 3). These storage time scales we determined agreed with 3D spherical spinel diffusion models we conducted at 1215°C for different crystal radii to assess the time scales of diffusive equilibration for spinels of different sizes (Fig. 1). The Mg compositions of Borgarhraun plagioclase macrocryst cores have equilibrated with the mush conditions (figs. S20 to S31). This means the coarsest of these crystals (1200-μm radius) required ~570 years of storage, assuming 3D diffusion in a spherical crystal at 1215°C and a core composition of An90 (fig. S32). The time scale of Mg in plagioclase equilibration provides a minimum estimate of the storage time, consistent with our spinel models that require hundreds to thousands of years of storage.

Fig. 3 Storage time estimates of near-Moho magma.

Cumulative frequency distributions showing the magma storage time estimates from the Bayesian inversion. Gray lines are distributions obtained by modeling individual spinel crystals. The black line is the cumulative frequency distribution of all of the modeled spinel crystals. The dashed blue line and the blue arrow are the minimum and the range of possible residence times based on Mg equilibration in plagioclase at 1215°C, respectively. The light-green region represents mush residence times based on bulk 3D spherical diffusion models shown in Fig. 1.

238U-230Th and 226Ra-230Th crystal residence ages from diverse tectonic settings range from 102 to 105 years (Fig. 4) (8). Discrepancies between these estimates and thermally activated diffusion time scales show that magmatic crystals may have complex and protracted histories throughout the crust that are largely independent of the melts that formed them and carried them to the surface for eruption (8). This highlights the difficulty in unpicking individual segments of magmatic evolution from bulk crystal ages integrated over zoned crystal populations. Overall, crystal residence time scales in basaltic systems have previously been estimated as hundreds to thousands of years by combining oxygen diffusion durations with U-series observations (27, 28), but these measurements do not provide a spatial context for the depth of magma residence. In contrast, textural and petrological constraints from Borgarhraun unambiguously link the diffusion time scale with near-Moho storage. Olivine-hosted melt inclusions link the crystal cargo to the carrier liquid by concurrent mixing and crystallization (11), meaning our results can be interpreted as magmatic storage times. Our diffusion chronometry approach therefore provides an opportunity for relative high temporal and spatial resolution of magmatic storage times in mafic lower-crustal systems (Fig. 4).

Fig. 4 Crystal residence times and storage depths for magmatic systems from different tectonic settings.

Shaded boxes represent the range of storage depth (normalized to the Moho) and residence time scale estimates made for individual volcanoes or volcanic systems. Time scale estimates for all eruptions except for Borgarhraun used U-series disequilibria, with residence times being calculated by using the difference between crystallization age and eruption age. The color of the box edges corresponds to the average SiO2 content of the whole rock sample used in the time scale estimates. The color of the symbols corresponds to the tectonic setting. Data and references are in table S4. kyr, thousand years.

Our storage time estimate for basaltic mushes in the lower crust have implications for understanding the time scale of magma transfer through the entire crustal column at a range of tectonic settings. This includes not only mid-ocean ridges and ocean islands (e.g., Iceland and Hawaii) but also oceanic and continental arc settings, where the lower crust is dominated by mafic and ultramafic magmatism, as shown by exposed crustal sections (e.g., the Kohistan arc in Pakistan) (29), primitive erupted products (5, 30), and geophysical observations (31). There are no direct constraints on the magmatic storage times in the lower crust of arcs. We therefore propose that our results provide an initial guide to expected near-Moho storage times in other mafic systems. For ocean island settings, our results indicate that primary magma can be stored near the Moho for thousands of years after initial mixing and crystallization (17). These magmas are rapidly transported over days to the mid- or shallow crust (7), where there is further long-term storage and magma processing (27, 28). This impacts our understanding of mantle melt transport, which is based on 226Ra excesses, and may lead to downward revision of melting region porosities when millennial storage of magma in the crust is considered (32).

Crystal storage in shallow, silicic arc systems can take place over time scales of up to 105 years, which is likely related to different thermal (1) and rheological properties when compared with basaltic systems. If the shorter processing times we have found apply to arc systems, then the upper crust may be episodically supplied with relatively young magma batches from the lower mafic crust (Fig. 4). The shorter crystal and magma storage times in the lower crust might reconcile differences between melt production rates in the mantle and the frequency of large silicic eruptions (1). In each tectonic setting, punctuated magma storage and rapid transport is likely a prominent feature. Although packets of magma fresh from the mantle melting region could take thousands of years to traverse the crust, storage could accommodate most of this time. This style of magma transfer will play a major role in the movement of heat and volatiles through the crust and controlling the connection between deep and surface geochemical reservoirs, such as the mantle and atmosphere.

Supplementary Materials

Materials and Methods

Supplementary Text

Figs. S1 to S34

Tables S1 to S4

References (3395)

References and Notes

Acknowledgments: We thank C. Richardson for helpful advice about FEniCS and M. Edmonds and D. Neave for their comments on an early draft of this manuscript. Funding: This research was funded by a NERC studentship awarded to E.J.F.M. (NE/L002507/1). Author contributions: The project was initially conceived by E.J.F.M. and J.M. E.J.F.M. developed the diffusion model. E.J.F.M. performed the EPMA and SEM analysis with help from I.B. I.B. conducted the QEMSCAN analysis and processing of QEMSCAN maps. T.J.B.H. conducted the thermodynamic modeling using THERMOCALC. E.J.F.M. did the diffusion modeling. E.J.F.M. and J.M. interpreted the results. E.J.F.M. wrote the first draft of the manuscript, which was revised by all authors. Competing interests: The authors declare that they have no competing financial interests. Data and materials availability: All data collected and modeled as part of this study are included in the supplementary materials and supplementary data files. The diffusion model code and Bayesian inversion in their current form are available from

Stay Connected to Science

Navigate This Article