Magmatic Densities Control Erupted Volumes in Icelandic Volcanic Systems

Magmatic density and viscosity exert fundamental controls on the eruptibility of magmas. In this study, we investigate the extent to which magmatic physical properties control the eruptibility of magmas from Iceland's Northern Volcanic Zone (NVZ). By studying subaerial flows of known age and volume, we are able to directly relate erupted volumes to magmatic physical properties, a task that has been near-impossible when dealing with submarine samples dredged from mid-ocean ridges. We find a strong correlation between magmatic density and observed erupted volumes on the NVZ. Over 85% of the total volume of erupted material lies close to a density and viscosity minimum that corresponds to the composition of basalts at the arrival of plagioclase on the liquidus. These magmas are buoyant with respect to the Icelandic upper crust. However, a number of small-volume eruptions with densities greater than typical Icelandic upper crust are also found in Iceland's neovolcanic zones. We use a simple numerical model to demonstrate that the eruption of magmas with higher densities and viscosities is facilitated by the generation of overpressure in magma chambers in the lower crust and uppermost mantle. This conclusion is in agreement with petrological constraints on the depths of crystallization under Iceland.

1. Introduction

It has long been understood that the range of lava compositions sampled at mid-ocean ridges is strongly linked to the buoyancy of magma in the oceanic crust (Huppert and Sparks, 1980; Sparks et al., 1980) and the viscosity of the magma (Walker, 1971). Primary melts generated at depth are expected to rise through the crust until ρ_m–ρ_c≈ 0, that is, the point at which the magma density ρ_m and crustal density ρ_c are equal. The ponding of magmas at their level of neutral buoyancy in the crust (Ryan, 1993) is considered to be an important factor in the prolonged existence of shallow magma chambers (Ryan, 1987).

Neutral buoyancy concepts can also be applied within magma chambers. Huppert and Sparks (1980) and Sparks et al. (1980) considered a fluid-dynamical approach to explore the conditions under which magma stored in stratified chambers may be erupted at the surface. As melt crystallizes within the chamber, the density of the residual liquid evolves until a density minimum is reached. Typical suites of mid-ocean ridge basalts (MORB) reach this density minimum at ~7–10 wt.% MgO (Sparks and Huppert, 1984); the buoyant residual liquids may then ascend to their level of neutral buoyancy, resulting in an eruption if the liquid remains buoyant all the way to the surface. Further crystallization of the residual liquid causes the density to increase once more, reducing the buoyancy and hence the eruptibility of the melt. Stolper and Walker (1980) defined a “window of eruptibility” for potential eruptive magmas, whereby the crust acts as an effective filter to high-density magmas. The range of compositions that can pass through this window depends firstly on the crustal density, and secondly on the position of the density minimum, which will itself vary depending on the initial composition and liquid line of descent of the parental melt. For mid-ocean ridge tholeiites erupting through typical oceanic crust, the window of eruptibility is found at melt MgO contents between 7 and 10 wt.%, or magnesium numbers (Mg#, molar Mg/(Mg+Fe²⁺) around 0.55–0.70. This fluid-dynamical approach predicts that basalts with MgO contents >10 wt.% are very unlikely to erupt (Huppert and Sparks, 1980). However, dense olivine-phyric magmas with whole-rock MgO contents >15 wt.% and Mg#>0.8 have been erupted on Iceland's neovolcanic zones. The eruption of these dense magmas cannot be simply explained by neutral buoyancy alone.

In this study, we explore the relationship between physical properties and erupted volumes for magmas from Iceland's Northern Volcanic Zone (NVZ). We then use a simple numerical model to explore how variations in magma chamber depth, country rock density and fracture strength, can facilitate the eruption of dense magmas to the surface.

2. Magmatism in Iceland

Iceland is situated where a mantle plume underlies the Mid-Atlantic Ridge. High mantle potential temperatures beneath Iceland (e.g., McKenzie and O'Nions, 1991; Langmuir et al., 1992; Matthews et al., 2016) cause a greater degree of melting than under normal mid-ocean ridges, producing anomalously thick crust that is exposed above sea level. Iceland is thus the only part of the global spreading system where lava compositions can be routinely tied to eruptions of known age and volume. Iceland's active neovolcanic zones provide an excellent opportunity to better understand the role of the crust as a density filter for magma. Previous studies of the density filter have focused on MORB, using statistical analysis of many hundreds of samples to determine the eruption probability for different magma compositions and densities. However, this approach is limited because it has rarely been possible to associate samples dredged from submerged spreading centers with individual eruptions of known age and volume. The geological controls available on Iceland allow us to study not only the range of erupted magma compositions but also the volumetric distribution of lavas with different compositions and different physical properties.

Neutral buoyancy and fluid-dynamical models assume that magma ascent will always be limited by neutral buoyancy. Eruption at the surface must therefore be facilitated by crystallization in one or more crustal magma chambers until the melt is sufficiently buoyant to be erupted at the surface, which is likely to occur when the magma is at or near the density minimum (Sparks et al., 1980). These models are appropriate for some volcanic settings in Iceland, notably central volcanoes such as Krafla where seismic surveys indicate the presence of a shallow magma chamber at ~3 km depth (Brandsdóttir et al., 1997). However, evidence from ophiolite sections (Kelemen et al., 1997; Korenaga and Kelemen, 1997) and from clinopyroxene crystallization pressures in basaltic lavas across Iceland (e.g., Maclennan et al., 2001) indicate that polybaric fractional crystallization occurs in sill-like bodies at a range of crustal depths, including near-Moho depths of 25–30 km (0.8–1.1 GPa; Maclennan et al., 2003; Maclennan, 2008; Winpenny and Maclennan, 2011; Neave and Putirka, 2017). The presence of modally layered gabbroic sills intruding dunite in the Moho transition zone of the Semail ophiolite suggests that magmas can be stored at, and undergo crystallization at, depths far greater than their level of magmatic neutral buoyancy (Kelemen et al., 1997) before the residual liquids ascend and erupt. Similarly, the presence of chrome-diopside in many of Iceland's lava shields suggests that these magmas did not undergo significant low-pressure fractionation, but were fed directly by primitive melts sourced from the base of the crust (Sigurdsson and Sparks, 1978). The lava shields include large-volume eruptions of more evolved basalts (e.g., Stóravíti, ~30 km³; Skjaldbreiður, ~11 km³) and smaller volumes of dense picrite containing up to 20% olivine macrocrysts (e.g., Háleyjabunga, 0.013 km³).

3. Magma Density and Viscosity

A dataset comprising 108 whole-rock samples from the Theistareykir and Krafla volcanic systems of Iceland's Northern Volcanic Zone was used to investigate the links between magmatic physical properties and the volumetric distribution of lava as a function of its composition. The samples correspond to postglacial eruptions of known volume and major element chemistry, and for each sample the modal proportion of macrocrysts has been visually estimated (Nicholson, 1990; Nicholson et al., 1991; Slater, 1996; Slater et al., 1998). For selected samples, we verified the modal proportions of olivine, plagioclase and rare clinopyroxene macrocrysts by point counting of thin sections. Macrocryst compositions in these sections were then determined using a Cameca SX100 electron microprobe at the University of Cambridge, using a beam diameter of 2 μm, current of 10 nA and accelerating voltage of 15 kV for all analyses. Whole-rock, point counting and mineral compositional data are provided as an electronic Datasheet 1.

The major element compositions of the NVZ samples are shown in Figure 1A. The change in trend at ~9.5 wt.% MgO reflects a change in fractionating mineral assemblage. Maclennan et al. (2001) used principal component analysis to demonstrate that the compositions of Theistareykir samples with 9.5–14.0 wt.% MgO can be explained by the addition or removal of a wehrlitic assemblage of 43% olivine, 47% clinopyroxene and 10% plagioclase, while the Theistareykir and Krafla samples with 6–9.5 wt.% MgO are best modeled by crystallization of a gabbroic assemblage with a modal mineralogy of ~10% olivine, ~40% clinopyroxene and ~50% plagioclase.

FIGURE 1

Figure 1. (A) Whole-rock major element compositions of postglacial lava samples from the Theistareykir (black circles) and Krafla (red diamonds) volcanic systems. Large filled symbols show samples that are aphyric, or for which the modal macrocryst content is not recorded. Large open symbols show samples with ≥1 vol.% macrocrysts. Small filled symbols show whole-rock compositions of Theistareykir and Krafla lavas of glacial or unknown age. (B) Relationship between whole-rock composition and density for postglacial NVZ samples, calculated assuming that the whole-rock composition represents a crystal-free magmatic liquid. Pale filled symbols with blue outlines show magmatic densities for macrocryst-bearing samples, recalculated from the liquid and mineral densities and volume fractions according to Equation (2).

The density of a melt consisting of N components with mole fractions x_j (j = 1,2,3…N) is typically calculated:

$\begin{array}{rc}\rho =\frac{{\displaystyle {\sum }_{j=1}^N}{x}_j{M}_j}{{\displaystyle {\sum }_{j=1}^N}{x}_j{V}_j}& (1)\end{array}$

where M is the molar mass and V the partial molar volume. This equation is often used to calculate the density of a magma from its whole-rock major element oxide composition, which is appropriate if the whole-rock composition represents a melt, i.e., a magmatic liquid. However, whole-rock analyses of phyric lavas provide the bulk average composition of the carrier melt and crystal components. Thus, the density calculated assuming that the whole-rock composition represents a magmatic liquid is not necessarily identical to the density of the magma, i.e., liquid plus crystals. For example, a melt carrying olivine macrocrysts will be denser than a liquid with the same bulk composition. Similarly, the liquid component of a magma might be buoyant with respect to the crust, but the same liquid carrying a cargo of olivine or clinopyroxene macrocrysts will have a higher bulk density and will not be able to ascend unless the macrocrysts settle from the liquid (Stolper and Walker, 1980). The density of a crystal-rich magma is best calculated as follows:

$\begin{array}{rc}{\rho }_m=F_{ol}{\rho }_{ol}+F_{pl}{\rho }_{pl}+(1-F_{ol}-F_{pl}){\rho }_l& (2)\end{array}$

where F is a volume fraction, ρ is a density, and the subscripts ol, pl, l, and m refer to olivine, plagioclase, liquid, and magma (i.e., liquid plus macrocrysts) respectively. Other macrocryst phases such as clinopyroxene may also be included in this calculation; however, clinopyroxene macrocrysts are very rare in the NVZ samples and so are not considered here.

We first calculated the density and viscosity of the NVZ samples assuming that the whole-rock composition represented a melt. We used Petrolog3 (Danyushevsky and Plechov, 2011) to calculate the liquidus temperature of the putative melt composition. The melt density at the liquidus temperature was then calculated from the molar volumes and mass fractions of the major element oxide components (Figure 1B). Partial molar volumes and thermal expansion coefficients were taken from Lange (1997), Lange and Carmichael (1990), Bottinga and Weill (1970), and Bottinga et al. (1982). Viscosities for the putative melt compositions were calculated at their liquidus temperatures using the model of Giordano et al. (2008).

To obtain a more realistic estimate of carrier liquid compositions in the NVZ, each whole-rock composition was corrected for the presence of macrocrysts. Point-counting data were used to obtain modal abundances for each macrocryst phase; where point-counting data were not available, we used published visual estimates of modal macrocryst abundance. Electron microprobe data were used to obtain representative compositions of each macrocryst phase in each sample. In individual samples, olivine and plagioclase macrocryst compositions tend to cluster around well-defined peaks with standard deviations of no more than 2.5 mol.% forsterite for olivine, and no more than 4 mol.% anorthite for plagioclase. Similar composition distributions have been obtained for olivine-phyric lavas across Iceland (Thomson and Maclennan, 2013). For each sample, the most probable macrocryst composition was used to calculate the mean macrocryst density, modeling forsterite-fayalite and albite-anorthite as ideal solid solutions with linear relationships between the densities of their endmembers (Fei, 1995). The mean macrocryst compositions were then subtracted from the whole-rock composition according to their modal abundances, and the remaining oxides renormalized to 100 wt.% to provide an estimate of the carrier liquid composition at the time of eruption. Densities and viscosities of these carrier liquids were calculated at their liquidus temperatures following the same method as for the whole-rock composition. Finally, the magmatic density (i.e., liquid plus crystals) was calculated from the liquid and mineral densities and volume fractions, according to Equation (2) (Figure 1B).

Of the 108 samples considered in this study, 63 were sparsely phyric to aphyric. Forty-five samples were recorded as containing >1% macrocrysts. Of these, 16 contained <5% macrocrysts and, for these samples, the difference between the density calculated from the whole-rock composition and the magmatic density (i.e., liquid plus crystals) is negligible (6 kg m⁻³ mean). The remaining samples contain up to 28.5% macrocrysts. These include olivine-phyric samples from Borgarhraun and the Theistareykir picrites, and plagioclase-phyric samples from Stóravíti. For the six samples containing ≥10 vol.% olivine macrocrysts, magmatic densities are on average 67 kg m⁻³ higher than the densities calculated under the assumption that their whole-rock compositions represent erupted liquids. For the ten samples containing ≥5 vol.% plagioclase macrocrysts and ≤2 vol.% of any other mineral phase, magmatic densities are on average 14 kg m⁻³ higher than the equivalent liquid density. The effect of uncertainty in the macrocryst modal proportions, either from point counting or visual estimation, on the calculated magmatic density is dependent on the density of the macrocryst phase of interest. An uncertainty of ±1% in the plagioclase modal abundance affects the calculated magmatic density by no more than ±1 kg m⁻³, while ±1% olivine affects the magmatic density by ±6 kg m⁻³.

Magmatic viscosities were calculated at the liquidus temperature of the carrier liquid from the melt viscosity and the macrocryst content, according to the Einstein-Roscoe equation:

$\begin{array}{rc}\mu ={\mu }_0{(1-\Phi /{\Phi }_m)}^{-n}& (3)\end{array}$

where μ₀ is the viscosity of the liquid and Φ is the crystal fraction. Φ_m and n are adjustable parameters that vary with the size, shape and distribution of particles in the melt; for magmatic processes these can be modeled using Marsh's constants Φ_m = 0.6 and n = 2.5 (Lejeune and Richet, 1995). We also tested the more sophisticated viscosity model of Costa (2005), which reduces to the Einstein-Roscoe equation at low crystal fractions. The two models yielded near-identical calculated magma viscosities for the NVZ samples, which is unsurprising given their low macrocryst contents.

At atmospheric pressure, calculated magmatic densities are at a minimum for samples with 0.635<Mg#<0.670. The density minimum roughly corresponds to the predicted arrival of plagioclase on the liquidus and the observed shift from dominantly olivine-phyric to plagioclase-phyric samples. The position of the density minimum is little affected by the choice of pressure at which the melt liquidus temperature is calculated: melt densities at the liquidus temperature of a given sample decrease by ~5–10 kg m⁻³ on increasing pressure from 0.001 to 0.9 GPa, and this effect is almost uniform for all of the samples considered. In order to highlight the control on physical properties from melt composition and macrocryst content, the results presented hereafter refer to densities calculated at atmospheric pressure.

4. Crystallization, Accumulation, and Vesiculation

The compositional array of magmas reflects the control of magmatic evolution by both fractional crystallization and melt mixing. The minimum available density and viscosity of magma is controlled by the fractionation and mixing paths of those magmas. The variation of density and viscosity as a function of melt fraction crystallized was calculated using Petrolog3. Calculations were run at a range of crustal pressures (0.001–0.9 GPa) and putative primary melt compositions. The choice of crystallization pressure affects the order in which olivine, plagioclase and clinopyroxene appear on the liquidus, but the position of the density minimum varies by only ~5–6 kg m⁻³. The density minimum is predicted to be lowest for isobaric crystallization at 0.9 GPa and highest for isobaric crystallization at 0.001 GPa, with polybaric crystallization models predicting density minima lying between the two isobaric cases. In all crystallization models we calculate a cusp in the predicted evolution of density and viscosity during fractional crystallization (Figure 2A). This cusp relates to the arrival of plagioclase on the basalt liquidus.

FIGURE 2

Figure 2. (A) Fractional crystallization of a NVZ picrite with a whole-rock composition in equilibrium with Fo₉₀ olivine. The crystallization model shown was calculated in Petrolog3 using a pressure of 0.4 GPa. The black solid line shows the density and viscosity evolution of the melt (i.e., residual liquid) assuming that it undergoes perfect fractional crystallization. The arrivals of plagioclase (P) and magnetite (M) on the liquidus are associated with shifts in the density evolution of the melt. Thin blue dashed lines show crystallization models for the same starting composition but under hydrous conditions, with the starting wt.% H₂O labeled on each line. The thick green dashed line shows the effect of adding olivine crystals to the starting composition on the magmatic (i.e., liquid plus crystals) density and viscosity; crosses show intervals of 10% olivine addition. (B) Relationship between magmatic density, viscosity, and erupted volume for postglacial lavas from Iceland's NVZ. The cumulative lava volume in each box is colored according to the scale at the top and contoured in 2 km³ intervals.

Subglacial basalts from the NVZ contain nominally undegassed H₂O contents of 0.1–0.5 wt.% (Nichols et al., 2002). The effect of adding up to 1 wt.% H₂O on melt liquid lines of descent was calculated in Petrolog3 using the mineral-melt equilibria of Danyushevsky (2001), and assuming that H₂O does not degas. Density and viscosity were calculated as a function of melt fraction crystallized, with water included as a phase in the density calculation using data from Ochs and Lange (1999). Increasing the H₂O content of the melt decreases its viscosity, and decreases the temperature at which olivine and plagioclase appear on the liquidus. The net result of these competing effects is that there is little change in the position of the viscosity minimum, where plagioclase joins the crystallizing assemblage (Figure 2A). Increasing the melt H₂O content decreases the melt density at all points along the crystallization pathway, but the position of the density minimum itself is minimally affected: for the NVZ picrite shown in Figure 2A the density minimum lies at 2,703 kg m⁻³ for an anhydrous melt, and 2,699 kg m⁻³ for a crystallization model with an initial 1 wt.% H₂O.

Volatile outgassing and vesiculation are expected to significantly lower the magmatic density and viscosity. Crystal-hosted melt inclusions from Borgarhraun indicate that nominally undegassed primitive magmas from the NVZ may contain up to 0.2 wt.% H₂O, 1,300 ppm CO₂ and 1,000 ppm S (Hauri et al., 2018). A basaltic melt with these volatile concentrations is predicted to become vapor-saturated at ~220 MPa (Newman and Lowenstern, 2002). Much of the CO₂ is expected to exsolve at pressures >50 MPa, and may be lost through open-system degassing if magma is stored in the shallow crust prior to eruption. After vapor saturation is reached H₂O exsolves continuously, but significant H₂O loss is not expected to occur until pressures around 10 MPa, while S may remain dissolved in the melt until pressures around 5 MPa. The pressure at which vigorous exsolution and bubble growth occurs is dependent on both the initial volatile content and the decompression rate, but is expected to occur once the melt has ascended to pressures <3 MPa (e.g., Gerlach, 1986). While the mean vesicularity in the NVZ samples considered in this study is 5% (range 0–25%), vesiculation is expected to occur primarily during magma ascent through the uppermost crust, and thus has minimal effect on magma buoyancy in the mid- to lower crust. Given that the volatile contents and decompression rates of the NVZ magmas in this study are poorly constrained, we hereafter assume volatile-free conditions.

5. Eruptive Volume Control by Physical Properties

The physical properties of magmas exert a strong control on the erupted volumes. By dividing the number of samples from a given eruption by the total volume of that eruption, we calculate a representative erupted volume for each NVZ sample. Hence, we obtain the volume of magma erupted at a given density and/or viscosity.

Despite the large range in predicted magma densities (2,670–2,840 kg m⁻³) and viscosities (5–180 Pa s) for postglacial lavas from the NVZ, the largest eruptions cluster close to the density and viscosity minima associated with the arrival of plagioclase on the basalt liquidus (Figure 2). Over 50% of the total erupted volume lies at densities <2,710 kg m⁻³, and over 85% below <2,720 kg m⁻³. Similarly, over 80% of the total erupted volume lies at viscosities <30 Pa s. The lowest density cluster on Figure 2B represents a single macrocryst-free sample from the postglacial lava shield Stóravíti (30 km³). Of the 27 samples from Stóravíti considered in this study, all but one have densities below 2,720 kg m³ and lie within the density well on Figure 2B. Stóravíti dominates the postglacial volumetric output of the Theistareykir volcanic system, but its exceptionally large volume cannot be simply explained by enhanced melt production and eruption rates during early postglacial times following deglaciation (Jull and McKenzie, 1996; Maclennan et al., 2002), since these are common to all the early postglacial eruptions in this study. Instead, we argue that it is the physical properties of the Stóravíti magma that result in the large volume of this eruption in comparison with other early postglacial eruptions such as the Theistareykir picrites.

Only 0.75% of the erupted volume from the NVZ has predicted density >2,750 kg m⁻³. This includes the Theistareykir picrites, and olivine-rich samples from Borgarhraun. Just over 5% of the total erupted volume has predicted viscosity >50 Pa s. Most of this volume corresponds to aphyric samples with the lowest MgO contents (MgO <7 wt.%) in the NVZ dataset; these higher viscosities therefore correspond to the samples with the lowest calculated liquidus temperatures rather than a population of magmas that have erupted after magnetite saturation and SiO₂ enrichment. The remaining high-viscosity samples (1% of the total erupted volume) have relatively low liquidus temperatures, and carry 6–16 vol.% macrocrysts. The volume maximum sits close to the minimum density and viscosity of magma available in the NVZ, while no samples in our dataset have both high density and high viscosity. This demonstrates the fundamental control of density and viscosity on volumetric output.

While upper crustal density plays an important role in controlling eruptive volumes, dense magmas (>2,800 kg m⁻³) are able to rise above their levels of neutral buoyancy in the crust. The eruption of such magmas is driven by the development of sufficient overpressure in near-Moho to mid-crustal magma chambers. One way to generate this overpressure is as buoyant melts encounter local permeability barriers such as unfractured or more rigid (including colder and thus stronger) country rock. The rapid crystallization of basaltic melt lenses also creates zones of reduced porosity and permeability such that ascending melts are trapped beneath these zones. Kelemen and Aharonov (1998) develop this model based on field evidence from the Moho transition zone of the Semail ophiolite, where gabbroic sills, relics of basaltic melt lenses, formed in denser harzburgitic country rock. Continuous melt influx to the region beneath a permability barrier leads to increasing overpressure, eventually generating a melt-filled fracture which may reach the surface and result in eruption. Alternatively, ascending magma may be halted by either encountering another permeability barrier, rheological contrast and/or rigidity contrast (e.g., Kavanagh et al., 2006; Menand, 2011), or by reaching its level of neutral buoyancy in the crust.

6. Linking Erupted Volumes to Physical Properties

6.1. Constant Overpressure in a Magma Chamber

The observed relationship between density, viscosity and erupted volume can be explored using simple models of the development of overpressure in melt lenses, which relate magma viscosity, density and input flux to ascent rate, volumetric output rate and eruption frequency (Kelemen and Aharonov, 1998). Melt influx to melt lenses at depth h may cause a fracture of width a to form when the overpressure exceeds the fracture stress of the country rock (Figure 3; see Table 1 for the notation used in the modeling). It is assumed that such fractures open instantaneously. Once open, the surrounding country rock is assumed to behave elastically with respect to the fracture. Melt flow through the fracture is assumed to be laminar and constant with height, and thus can be described by an average melt velocity. Assuming constant overpressure and ignoring horizontal extensional stresses, the dyke width and melt flow speed are given by:

$\begin{array}{rc}a=\frac{h{P}_f}{G},\omega \approx \frac{a^2}{12\mu h}(P_f+\Delta \rho gh)& (4)\end{array}$

The predicted flux, the product aω, was calculated using the above equations for a range of magma densities and viscosities, setting P_f = 1 MPa, G = 10 GPa and h = 10 km (Figure 4). The magma chamber depth is set to lie within the Icelandic mid-crust. The values of P_f and G are not well known. Kelemen and Aharonov (1998) selected a value of G = 10 GPa after Turcotte and Schubert (2014), although Auriac et al. (2013) suggest that G for the Icelandic lower crust may lie closer to 75 GPa. P_f, the overpressure required for fracture, is poorly constrained. Kelemen and Aharonov (1998) selected a value of 10 MPa, but we initally found a good match between our observations and model predictions at P_f = 1 MPa.

FIGURE 3

Figure 3. Schematic illustration of a crustal magma reservoir at depth h, with an open fracture allowing magma to ascend with velocity ω from the top of the reservoir. Modified after Kelemen and Aharonov (1998). The fracture of width a forms when the overpressure in the magma chamber exceeds the fracture stress of the crust. The magma reservoir may be fed directly by mantle melts or from deeper, sill-like melt lenses. Under certain conditions, high-density olivine-phyric magmas are enabled to erupt from deep magma chambers directly to the surface. See text for details.

TABLE 1

Table 1. Parameters used in eruption modeling.

FIGURE 4

Figure 4. Comparison of sill model predictions and observed relationship between erupted volume and magmatic density and viscosity. (A) Cumulative histogram of the relationship between inverse viscosity and fraction of total erupted volume. The black continuous line shows all observations from the NVZ. The blue line shows NVZ samples with density under 2,720 kg m⁻³. The predicted relationship between magma viscosity and maximum melt flux rate is shown as a dashed red curve. This curve was calculated only for magma with the viscosity range in samples with densities under 2,720 kg m⁻³. The density of the model magma was held constant at 2,710 kg m⁻³ for these calculations. (B) Cumulative histogram showing the relationship between erupted volume fraction and sample density, relative to a reference density of 2,710 kg m⁻³. The black continuous line shows all observations from the NVZ; the blue line shows samples with viscosity under 50 Pa s. The model results shown as a dashed red curve were calculated for the density range found in melts with viscosity under 50 Pa s. The viscosity of the model magma was held constant at 50 Pa s for these calculations. For both plots, the red curves are calculated using the following reference values: G = 10 GPa; h = 10 km; P_f = 1 MPa. The green dashed curve in (B) shows the effect of changing the reference values to G = 75 GPa and P_f = 3 MPa.

This generalized model cannot encapsulate all the details of individual eruptions (for example, the assumption of constant a is unlikely to be correct) but illustrates possible mechanisms for the control exerted by physical properties. Given its simplicity, it provides a surprisingly good fit to the observed relationship between physical properties and erupted volumes. There is some trade-off between P_f, G, and h as we attempt to fit the NVZ observations: for example, higher values of G will require higher values of P_f. Model predictions match the NVZ observations with a wide range of values for G and h: there is little difference in the model predictions for 10< G <100 GPa. The model predictions are most sensitive to variations in P_f. We obtain a reasonable fit to the observations for P_f values up to 3. For P_f > 3 the model predicts that dense magmas will be buoyant, but the model curves no longer fit the observations. We emphasize that, while the details of our model predictions are unlikely to be correct, our model nonetheless facilitates a simple physical exploration of the plausibility of the suggested relationship between erupted volumes and magmatic physical properties.

The influence of viscosity on the erupted volumes can be isolated when the magma density is held fixed in the models, or the observations filtered to lie within a narrow density range. Both model predictions are predominantly concave-up on Figure 4A, indicating that greater erupted volumes are found at lower viscosity. The role of buoyancy is clear from inspection of the results of calculations performed with fixed magma viscosity and varying density (Figure 4B). These model predictions indicate that as a result of the magmatic overpressure ~15% of the erupted melt may have a density greater than that of the country rock. Assuming a reasonable country rock density of 2,710 kg m⁻³ (e.g., Staples et al., 1997), this proportion is roughly in agreement with the observations. The model results predict that magmas with density >10–15 kg m⁻³ greater than that of the country rock are not erupted. However, 5–10% of the observed erupted volume is composed of magma with densities of 2,730 kg m⁻³ or greater, and some picrites have densities of >2,800 kg m⁻³, well in excess of expected Icelandic upper crustal densities. The eruption of such magma is not predicted by the simple sill models and may indicate that magma chambers feed eruptions from a variety of physical surroundings, as explored below.

6.2. Constant Volumetric Influx Rate to Magma Chamber

In order to account for the eruption of dense magmas through low-density upper crust, the physical conditions of the magma chambers supplying Icelandic eruptions must vary. More complex models were used to investigate the role of varying magma chamber depth, country rock density and fracture strength. This fuller realization of the models was explored by solving the differential equations for inflation of a magma chamber, fracture of its walls, release of magmatic overpressure and solidification of the dyke. Rather than maintaining constant magmatic overpressure, this scheme involves constant volumetric influx rate to the magma chamber and allows significant overpressure fluctuations to occur. The governing equations, posed by Kelemen and Aharonov (1998), were solved with a 4th order Runge-Kutta scheme paying particular attention to numerical accuracy at the time of fracture development, when predicted magma fluxes vary rapidly. The output magmatic flux was compared with an input at constant rate to the chamber.

Increasing h (i.e., deeper magma chamber), decreasing P_f, and/or increasing G encourage eruption at high densities and viscosities (Figure 5). Therefore, increased chamber depth may allow picrites to erupt at the surface, past their level of neutral buoyancy. Similarly, if the density of the country rock is increased, then eruption of denser magma is permitted. Magma chambers for picrites are therefore likely to be in the lower crust or uppermost mantle, where country rock densities of ≥2,950 kg m⁻³ are expected (Staples et al., 1997). This inference is in agreement with petrological constraints on the depth of magma chambers, which indicate that crystallization of picrites and high-Mg# basalts under Theistareykir occurs at depths of 15–30 km (Maclennan et al., 2001, 2003; Winpenny and Maclennan, 2011).

FIGURE 5

Figure 5. Results of constant volumetric influx rate models. Lines show boundary between density-viscosity corresponding to no predicted eruption (top right) and predicted eruption (bottom left). The solid black line shows an arbitrary reference set of conditions, with magma chamber depth at h = 20 km, P_f = 5 MPa, and G = 25 GPa. Note that eruption is expected for a wide range of magma densities and viscosities from this chamber. The other lines show the effect of varying the parameters individually, all in a sense which reduces the range of melt densities and viscosities that can be erupted. Dotted black line: effect of reducing P_f to 4 MPa, with all other parameters set at the reference values. Solid gray line: effect of setting h = 10 km with all other parameters set at the reference values. Dashed black line: effect of increasing G to 50 GPa, with all other parameters set at the reference values.

7. Conclusions

The physical properties of magmas exert key controls on the erupted volumes from Iceland's active neovolcanic zones. The largest volume of erupted material sits at a well-defined density and viscosity minimum coincident with the arrival of plagioclase on the basalt liquidus. These basalts are buoyant with respect to the Icelandic upper crust. Almost 70% of the erupted volume in the NVZ comprises basalts that sit close to this minimum.

The density filter is, however, leaky. Several small-volume eruptions with densities greater than the Icelandic upper crust are found in the NVZ. These eruptions include primitive basalts and olivine-phyric picrites, which may contain up to ~25% olivine macrocrysts. The eruption of these dense magmas is facilitated by the generation of overpressure in magma chambers in the lower crust or uppermost mantle. This conclusion is consistent with petrological constraints on the depth of crystallization of high-Mg# melts under Iceland.

Author Contributions

Initial concept by JM. MH conducted the point-counting and electron microprobe analyses; calculated magmatic densities, viscosities and crystallization models; and wrote the manuscript. JM provided samples and wrote the numerical code solving the Kelemen and Aharonov (1998) governing equations. Both authors contributed to data interpretation.

Funding

We acknowledge support from the Department of Earth Sciences, University of Cambridge, and funding through NERC grants NE/N018575/1 and NE/P002331/1.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Thierry Menand and Oli Shorttle for their comments on earlier versions of this work, and two reviewers for their thoughtful comments which improved and clarified the manuscript. We thank Olivier Bachmann for editorial handling.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2018.00029/full#supplementary-material

Datasheet 1. (1) Eruption volumes, whole-rock compositions, macrocryst abundances, and calculated densities for samples from Iceland's Northern Volcanic Zone. (2) Mineral point-counting data for NVZ postglacial basalt samples. (2) Mineral compositional data for NVZ postglacial basalts.

References

Auriac, A., Spaans, K. H., Sigmundsson, F., Hooper, A., Schmidt, P., and Lund, B. (2013). Iceland rising: solid Earth response to ice retreat inferred from satellite radar interferometry and visocelastic modeling. J. Geophys. Res. Solid Earth 118, 1331–1344. doi: 10.1002/jgrb.50082

CrossRef Full Text | Google Scholar

Bottinga, Y., Weill, D., and Richet, P. (1982). Density calculations for silicate liquids. I. Revised method for aluminosilicate compositions. Geochim. Cosmochim. Acta 46, 909–919. doi: 10.1016/0016-7037(82)90047-3

CrossRef Full Text | Google Scholar

Bottinga, Y., and Weill, D. F. (1970). Densities of liquid silicate systems calculated from partial molar volumes of oxide components. Amer. J. Sci. 269, 169–182. doi: 10.2475/ajs.269.2.169

CrossRef Full Text | Google Scholar

Brandsdóttir, B., Menke, W., Einarsson, P., White, R. S., and Staples, R. K. (1997). Faroe-Iceland ridge experiment. 2. Crustal structure of the Krafla central volcano. J. Geophys. Res. Solid Earth 102, 7867–7886. doi: 10.1029/96JB03799

CrossRef Full Text | Google Scholar

Costa, A. (2005). Viscosity of high crystal content melts: dependence on solid fraction. Geophys. Res. Lett. 32:L22308. doi: 10.1029/2005GL024303

CrossRef Full Text | Google Scholar

Danyushevsky, L. V. (2001). The effect of small amounts of H₂O on crystallisation of mid-ocean ridge and backarc basin magmas. J. Volcanol. Geother. Res. 110, 265–280. doi: 10.1016/S0377-0273(01)00213-X

CrossRef Full Text | Google Scholar

Danyushevsky, L. V., and Plechov, P. (2011). Petrolog3: integrated software for modeling crystallization processes. Geochem. Geophys. Geosyst. 12, 1–32. doi: 10.1029/2011GC003516

CrossRef Full Text | Google Scholar

Fei, Y. (1995). “Thermal expansion,” in Mineral Physics and Crystallography: A Handbook of Physical Constants, vol. 2, ed T. J. Ahrens (Washington, DC: AGU), 29–44.

Google Scholar

Gerlach, T. M. (1986). Exsolution of H₂O, CO₂, and S during eruptive episodes at Kilauea Volcano, Hawaii. J. Geophys. Res. 91, 12177–12185. doi: 10.1029/JB091iB12p12177

CrossRef Full Text | Google Scholar

Giordano, D., Russell, J. K., and Dingwell, D. B. (2008). Viscosity of magmatic liquids: a model. Earth Planet. Sci. Lett. 271, 123–134. doi: 10.1016/j.epsl.2008.03.038

CrossRef Full Text | Google Scholar

Hauri, E. H., Maclennan, J., McKenzie, D., Gronvold, K., Oskarsson, N., and Shimizu, N. (2018). CO₂ content beneath northern Iceland and the variability of mantle carbon. Geology 46, 55–58. doi: 10.1130/G39413.1

CrossRef Full Text | Google Scholar

Huppert, H. E., and Sparks, R. S. J. (1980). Restrictions on the compositions of mid-ocean ridge basalts - A fluid dynamical investigation. Nature 286, 46–48. doi: 10.1038/286046a0

CrossRef Full Text | Google Scholar

Jull, M., and McKenzie, D. (1996). The effect of deglaciation on mantle melting beneath Iceland. J. Geophys. Res. Solid Earth 101, 21815–21828. doi: 10.1029/96JB01308

CrossRef Full Text | Google Scholar

Kavanagh, J. L., Menand, T., and Sparks, R. S. J. (2006). An experimental investigation of sill formation and propagation in layered elastic media. Earth Planet. Sci. Lett. 245, 799–813. doi: 10.1016/j.epsl.2006.03.025

CrossRef Full Text | Google Scholar

Kelemen, P. B., and Aharonov, E. (1998). “Periodic formation of magma fractures and generation of layered gabbros in the lower crust beneath oceanic spreading ridges,” in Geophysical Monograph, vol. 106, eds W. R. Buck, P. T. Delaney, J. A. Karson, and Y. Lagabrielle (Washington, DC: AGU), 267–289.

Google Scholar

Kelemen, P. B., Koga, K., and Shimizu, N. (1997). Geochemistry of gabbro sills in the crust-mantle transition zone of the Oman ophiolite: implications for the origin of the oceanic lower crust. Earth Planet. Sci. Lett. 146, 475–488. doi: 10.1016/S0012-821X(96)00235-X

CrossRef Full Text | Google Scholar

Korenaga, J., and Kelemen, P. B. (1997). Origin of gabbro sills in the Moho transition zone of the Oman ophiolite: implications for magma transport in the oceanic lower crust. J. Geophys. Res. 102, 27729–27749. doi: 10.1029/97JB02604

CrossRef Full Text | Google Scholar

Lange, R. A. (1997). A revised model for the density and thermal expansivity of K₂O-Na₂O-CaO-MgO-Al₂O₃-SiO₂ liquids from 700 to 1900 K: extension to crustal magmatic temperatures. Contrib. Mineral. Petrol. 130, 1–11. doi: 10.1007/s004100050345

CrossRef Full Text | Google Scholar

Lange, R. L., and Carmichael, I. S. E. (1990). Thermodynamic properties of silicate liquids with emphasis on density, thermal-expansion and compressibility. Rev. Mineral. Geochem. 24, 25–64.

Google Scholar

Langmuir, C. H., Klein, E. M., and Plank, T. (1992). “Petrological constraints on melt formation and migration beneath mid-ocean ridges,” in Geophysical Monograph, vol. 71, eds J. Phipps Morgan, D. Blackman, and J. M. Sinton (American Geophysical Union), 183–280.

Lejeune, A. M., and Richet, P. (1995). Rheology of crystal-bearing silicate melts - an experimental study at high viscosities. J. Geophys. Res. Solid Earth 100, 4215–4229. doi: 10.1029/94JB02985

CrossRef Full Text | Google Scholar

Maclennan, J. (2008). Concurrent mixing and cooling of melts under Iceland. J. Petrol. 49, 1931–1953. doi: 10.1093/petrology/egn052

CrossRef Full Text | Google Scholar

Maclennan, J., Jull, M., McKenzie, D., Slater, L., and Grönvold, K. (2002). The link between volcanism and deglaciation in Iceland. Geochem. Geophys. Geosyst. 3:1062. doi: 10.1029/2001GC000282

CrossRef Full Text | Google Scholar

Maclennan, J., McKenzie, D., Grönvold, K., Shimizu, N., Eiler, J. M., and Kitchen, N. (2003). Melt mixing and crystallization under Theistareykir, northeast Iceland. Geochem. Geophys. Geosyst. 4:8624. doi: 10.1029/2003GC000558

CrossRef Full Text | Google Scholar

Maclennan, J., McKenzie, D., Grönvold, K., and Slater, L. (2001). Crustal accretion under northern Iceland. Earth Planet. Sci. Lett. 191, 295–310. doi: 10.1016/S0012-821X(01)00420-4

CrossRef Full Text | Google Scholar

Matthews, S., Shorttle, O., and Maclennan, J. (2016). The temperature of the Icelandic mantle from olivine-spinel aluminum exchange thermometry. Geochem. Geophys. Geosyst. 17, 4725–4752. doi: 10.1002/2016GC006497

CrossRef Full Text | Google Scholar

McKenzie, D., and O'Nions, R. K. (1991). Partial melt distributions from inversion of rare-earth element concentrations. J. Petrol. 32, 1021–1091. doi: 10.1093/petrology/32.5.1021

CrossRef Full Text | Google Scholar

Menand, T. (2011). Physical controls and depth of emplacement of igneous bodies: a review. Tectonophysics 500, 11–19. doi: 10.1016/j.tecto.2009.10.016

CrossRef Full Text | Google Scholar

Neave, D. A., and Putirka, K. D. (2017). A new clinopyroxene-liquid barometer, and implications for magma storage pressures under Icelandic rift zones. Amer. Mineral. 102, 777–794. doi: 10.2138/am-2017-5968

CrossRef Full Text | Google Scholar

Newman, S., and Lowenstern, J. B. (2002). VolatileCalc: a silicate melt-H₂O-CO₂ solution model written in Visual Basic for excel. Comput. Geosci. 28, 597–604. doi: 10.1016/S0098-3004(01)00081-4

CrossRef Full Text | Google Scholar

Nichols, A. R. L., Carroll, M. R., and Hóskuldsson, A. (2002). Is the Iceland hot spot also wet? Evidence from the water contents of undegassed submarine and subglacial pillow basalts. Earth Planet. Sci. Lett. 202, 77–87. doi: 10.1016/S0012-821X(02)00758-6

CrossRef Full Text | Google Scholar

Nicholson, H. (1990). The Magmatic Evolution of Krafla, NE Iceland. Ph.D. thesis, University of Edinburgh.

Nicholson, H., Condomines, M., Fitton, J. G., Fallick, A. E., Grönvold, K., and Rogers, G. (1991). Geochemical and isotopic evidence for crustal assimilation beneath Krafla, Iceland. J. Petrol. 32, 1005–1020. doi: 10.1093/petrology/32.5.1005

CrossRef Full Text | Google Scholar

Ochs, F. A., and Lange, R. A. (1999). The density of hydrous magmatic liquids. Science 283, 1314–1317. doi: 10.1126/science.283.5406.1314

PubMed Abstract | CrossRef Full Text | Google Scholar

Ryan, M. P. (1987). “Neutral buoyancy and the mechanical evolution of magmatic systems,” in Magmatic Processes: Physicochemical Principles, ed B. O. Mysen (Geochemical Society Special Publication 1), 259–287.

Google Scholar

Ryan, M. P. (1993). Neutral buoyancy and the structure of mid-ocean ridge magma reservoirs. J. Geophys. Res. Solid Earth 98, 22321–22338. doi: 10.1029/93JB02394

CrossRef Full Text | Google Scholar

Sigurdsson, H., and Sparks, R. S. J. (1978). Lateral magma flow within rifted Icelandic crust. Nature 274, 126–130. doi: 10.1038/274126a0

CrossRef Full Text | Google Scholar

Slater, L. (1996). Melt Generation Beneath Iceland. Ph.D. thesis, University of Cambridge.

Google Scholar

Slater, L., Jull, M., McKenzie, D., and Grönvold, K. (1998). Deglaciation effects on mantle melting under Iceland: results from the Northern Volcanic Zone. Earth Planet. Sci. Lett. 164, 151–164. doi: 10.1016/S0012-821X(98)00200-3

CrossRef Full Text | Google Scholar

Sparks, R. S. J., and Huppert, H. E. (1984). Density changes during the fractional crystallization of basaltic magmas - Fluid dynamic implications. Contrib. Mineral. Petrol. 85, 300–309. doi: 10.1007/BF00378108

CrossRef Full Text | Google Scholar

Sparks, R. S. J., Meyer, P., and Sigurdsson, H. (1980). Density variation amongst mid-ocean ridge basalts - Implications for magma mixing and the scarcity of primitive lavas. Earth Planet. Sci. Lett. 46, 419–430. doi: 10.1016/0012-821X(80)90055-2

CrossRef Full Text | Google Scholar

Staples, R. K., White, R. S., Brandsdóttir, B., Menke, W., Maguire, P. K. H., and McBride, J. H. (1997). Faroe-Iceland Ridge Experiment. 1. Crustal structure of northeastern Iceland. J. Geophys. Res. Solid Earth 102, 7849–7866. doi: 10.1029/96JB03911

CrossRef Full Text | Google Scholar

Stolper, E., and Walker, D. (1980). Melt density and the average composition of basalt. Contrib. Mineral. Petrol. 74, 7–12. doi: 10.1007/BF00375484

CrossRef Full Text | Google Scholar

Thomson, A., and Maclennan, J. (2013). The distribution of olivine compositions in Icelandic basalts and picrites. J. Petrol. 54, 745–768. doi: 10.1093/petrology/egs083

CrossRef Full Text | Google Scholar

Turcotte, D., and Schubert, G. (2014). Geodynamics, 3rd Edn. Cambridge University Press.

Walker, G. P. L. (1971). Viscosity control of the composition of ocean floor volcanics. Philos. Trans. R. Soc. Ser. A 268, 727–729. doi: 10.1098/rsta.1971.0023

CrossRef Full Text | Google Scholar

Winpenny, B., and Maclennan, J. (2011). A partial record of mixing of mantle melts preserved in icelandic phenocrysts. J. Petrol. 52, 1791–1812. doi: 10.1093/petrology/egr031

CrossRef Full Text | Google Scholar