Continuous mush disaggregation during the long-lasting Laki fissure eruption, Iceland

Macrocryst entrainment from disaggregating crystal mushes is a widely recognized process in basaltic plumbing systems (Rhodes et al. 1979; Lange et al. 2013a; Neave et al. 2014). For example, primitive olivine crystals carried by evolved liquids have been interpreted as entrained macrocrysts in samples from Iceland (Thomson and Maclennan 2013), Hawai'i (Vinet and Higgins 2010), Réunion (Albarède and Tamagnan 1988), and mid-ocean ridges (Donaldson and Brown 1977). Basalt-hosted high-anorthite plagioclase macrocrysts have also been understood as disaggre-gated mush remnants in a range of geological settings (Hansen and Grönvold 2000; Ridley et al. 2006). Moreover, isotope and trace element disequilibria between entrained macrocrysts and their carrier liquids indicate that crystals and melts are often derived from different mantle melt distributions (Halldórsson et al. 2008; Winpenny and Maclennan 2011; Lange et al. 2013b). Mush disaggregation has also been inferred from whole-rock geochemical systematics, whereby elemental abundances in variably porphyritic bulk samples are controlled by element compatibilities in mush-derived components (Salaün et al. 2010; Passmore et al. 2012). Finally, macrocryst entrainment has been identified in rock textures, from the size, shape, and abundance of macrocrysts themselves (Higgins 1996; Holness et al. 2007).

However, converting information from thin sections into CSDs faces two major challenges. First, to plot 2D crystal size data obtained from thin sections on explicitly 3D CSD plots it must be stereologically corrected. While Higgins (2000) presented an internally consistent tool for such corrections (CSDCorrections), 3D crystal shapes need to be uniform and well-defined to obtain reliable results, which is often not the case. Second, generating data sets for CSD calculations is time-consuming. Before crystal shape data can be determined using image analysis tools, thin section images must be digitized, which usually involves extensive manual input (Higgins 2000; Shea et al. 2010). Instruments capable of automated phase mapping (e.g., Gottlieb et al. 2000; Pirrie et al. 2004) thus represent an appealing, but as yet unevaluated, tool for the rapid acquisition of images suitable for CSD production.

In this contribution, we present phase proportion and plagioclase size data acquired from the products of the long-lasting 1783–1784 CE Laki eruption in Iceland using manual and automated Quantitative Evaluation of Minerals by SCANing electron microscopy (QEMSCAN) approaches. Using these data, we first evaluate whether automated approaches provide the same textural information as manual approaches in a reproducible and thus reliable way. We then investigate whether mush disaggregation occurred at the onset of the Laki eruption or throughout its eight-month duration, which has important implications for understanding how large volumes of magma are mobilized in basaltic fissure eruptions and interpreting future signs of unrest.

The 1783–1784 CE Laki eruption, also known as the Skaftáreldar (Skaftár Fires), took place over eight months between June 8, 1783 and February 7, 1784, in the Síða highlands of southern Iceland as part of a two-year volcano-tectonic episode in the Grímsvötn volcanic system (Fig. 1). In total, 15.1 km³ of lava and tephra erupted along a 27 km long series of 10 en echelon fissures that opened progressively from the southwest to the northeast forming Lakagígar crater row (Thordarson and Self 1993). The opening of many fissures was preceded by elevated levels of seismicity that marked the onset of new eruptive episodes (Thordarson and Self 1993). Each of the 10 eruptive episodes is thought to have started with a brief period of explosive activity that transitioned via fire fountaining to lava effusion (Guilbaud et al. 2007). Fissures I–V opened to the southwest of Laki Mountain and discharged primarily down the Skaftá river gorge, whereas fissures VI–X opened to the northeast of Laki Mountain and discharged primarily down the Hverfisfljót river gorge. The extensive lava field produced during the eruption covers ~600 km² and extends up to 40 km from the source vents. Thermally efficient transport over these long distances resulted in continued evolution of the lava during emplacement (Guilbaud et al. 2007).

Although Sigmarsson et al. (1991) reported extreme homogeneity in the isotope and trace element composition of the Laki lava flow, Passmore et al. (2012) subsequently noted subtle but statistically significant variations in the composition of whole-rock samples that correlate with their macrocryst contents. Passmore et al. (2012) thus proposed that whole-rock variability in the Laki flow reflects the presence of varying amounts of a mush-derived component containing incompatible element-poor macrocrysts and incompatible element-rich interstitial melts. Using a combination of petrography and microanalysis, Neave et al. (2013) subsequently identified an assemblage of primitive macrocryst cores that were interpreted as the crystal component of this disaggregated mush.

Four samples representing a range of eruptive episodes were selected from the 54 basaltic specimens described by Passmore et al. (2012) to investigate disaggregation processes throughout course of the eruption. Sampling locations are shown in Figure 1. With the exception of distally collected LAK27, which represents a crystal-rich end-member, we chose proximal samples (<1–15 km from the vent) to minimize the effects of textural evolution during transport. These samples were chosen to be representative of their episodes on the basis of their petrography and total crystal contents, which are within a few percent (absolute) of the mean values for their respective episodes (Passmore et al. 2012).

Two samples were selected from each of episodes I–V (LAK18 and LAK09) and episodes VI–X (LAK04 and LAK27). LAK18 and LAK09 were erupted southwest of Laki Mountain during episodes I and III, respectively, the highest mass flux episodes of the eruption (Thordarson and Self 1993). These rapidly quenched, porphyritic samples have low crystal contents (<10% by area) and were collected close to their source vents (~5 and <1 km, respectively; Figs. 1 and 2c; Supplementary¹Figs. 1a and 1b). LAK04 and LAK27 contain more crystals (>10% by area) and were erupted from fissures VII and VIII, respectively to the northeast of Laki Mountain (Fig. 1). LAK04 is porphyritic and was collected at a moderate distance from its source vent (~15 km) in the Skaftá river gorge where it cooled sufficiently rapidly to form a glassy rind (Fig. 2d; Supplementary¹Fig. 1c). In contrast, seriate LAK27 was collected from a distal flow lobe on the Síða Plain (~40 km) and has a coarse groundmass (Figs. 2a and 2b; Supplementary¹Fig. 1d).

All samples contain crystals of plagioclase, clinopyroxene, and olivine (Figs. 2a and 2b). To retain consistency with previous studies, the following apparent long-axis lengths (L) were used as thresholds for excluding groundmass grains: L = 0.15 mm for plagioclase and L = 0.20 mm for clinopyroxene and olivine (Passmore et al. 2012). Plagioclase crystals with L ≥ 0.15 mm were further divided into two populations: smaller microcrysts (0.15 > L < 0.5 mm) and larger macrocrysts (L > 0.5 mm). Within the suite of four samples investigated here, plagioclase macrocrysts reach up to ~3 mm in length, while clinopyroxene and olivine macrocrysts reach only ~1.5 mm. In other samples from Laki, plagioclase and clinopyroxenes macrocrysts can reach up to ~8 mm. Both micro- and macrocrysts have systematically more primitive compositions than groundmass grains (Neave et al. 2013): >An₆₅, Mg#_Cpx > 75, and >Fo₇₁.

Point-counting data used in this study were collected by Passmore et al. (2012), who also provided a detailed methodological description. In summary, each slide was point-counted 3–5 times at the School of GeoSciences, University of Edinburgh, U.K., using a manually operated mechanical slide holder moved in <0.20 mm increments in x and y directions. The following phases were counted: plagioclase, clinopyroxene, olivine, vesicles, and groundmass, which encompassed glass, mesostasis, and groundmass grains. Over 1000 points were counted in each repeat, and phase proportions were determined with the following 1σ relative precisions according to repeat measurements: ±18.4% for plagioclase, ±33.6% for clinopyroxene, ±45.5% for olivine, and ±6.0% for the total crystal content.

QEMSCAN images were produced using a Quanta 650F, field emission gun (FEG) scanning electron microscope (SEM), equipped with two Bruker XFlash 6130 energy-dispersive X-ray spectrometers (EDS) at the Department of Earth Sciences, University of Cambridge, U.K. The fully automated system includes an energy-dispersive X-ray (EDX) spectrum acquisition and classification procedure. Analyses were performed by obtaining field scans that provided a complete characterization of sample (thin section) surfaces above a predefined backscattered electron (BSE) threshold (Gottlieb et al. 2000; Pirrie et al. 2004). BSE brightness coefficients used to apply this threshold were calibrated against quartz, gold, and copper standards. EDX spectra for each pixel were generated from 2000 X-ray counts at 25 kV and 10 nA, and at spatial resolutions (i.e., pixel sizes) of 200, 50, and 10 μm. Imaging times ranged from 2 min for a resolution of 200 μm to 8 h for a resolution of 10 μm. Spectra were then processed using species identification protocol (SIP) files that discriminated user-defined minerals on the basis of their characteristic X-ray and electron backscatter intensities computed from ideal mineral compositions normalized to the beam conditions. Images collected at spatial resolutions of 200, 50, and 10 μm are henceforth referred to as Q200, Q50, and Q10 images, respectively. Q10 images are shown in Figure 3, and full images collected at all resolutions are provided in the electronic appendix¹.

SIP files were produced using quantitative electron microprobe analyses from the same samples for guidance (Passmore et al. 2012; Neave et al. 2013). Plagioclase zoning was resolved using SIP files that sorted plagioclase pixels into the following compositional bins based on their relative Si, Al, Ca, and Na contents: anorthite (~An₉₀), bytownite (~An₈₀), labradorite (~An₆₀), and andesine (~An₄₀). To improve the textural classification of plagioclase grains, images were also processed with SIP files tuned to distinguish groundmass pixels (<An₆₅) from micro- and macrocryst pixels (>An₆₅). Comparing QEMSCAN images with published electron microprobe data suggests that plagioclase anorthite contents were generally determined with a precision better than 20 mol% under the conditions used; QEMSCAN images nevertheless contain some noise reflecting the incorrect binning of some pixels. Olivine zoning was resolved in a similar manner by using relative Mg and Fe contents.

Plagioclase sizes were obtained from two suites of digitized images. The first suite of images was obtained by manually tracing plagioclase crystals on high-resolution (4000 dpi) thin sections scans in Inkscape, which took approximately 20 h per thin section. For consistency with point-counting data, we excluded ground-mass grains with L < 0.15 mm. Vesicles were also traced to calculate vesicle-free population densities. Crystal and vesicle size data then were extracted from these digitized images using the analyze particles tool in ImageJ (Abràmoff et al. 2004).

The second suite of images was generated from QEMSCAN analyses using FEI's iExplorer software package. Plagioclase sizes were estimated from these images in two ways. First, plagioclase sizes were determined using the granulator tool within iExplorer that segments and measures discrete particles of user-specified phases. Second, plagioclase sizes were determined from a set of granulated QEMSCAN images that were manually rectified to separate glomerocrysts into individual macrocrysts (e.g., Shea et al. 2010). Rectification was necessary because QEMSCAN analyses do not discriminate between touching crystals of the same composition; QEMSCAN images are phase maps, not grain maps. Rectified images were then measured using ImageJ in the same manner as the manually traced images.

Plagioclase areas determined using ImageJ and iExplorer were calculated by summing the total number of pixels within each particle in images to which thresholds had been applied. Uncertainties in particle areas therefore reflect a trade-off between image resolution and particle size. Given that each erroneous pixel in an 80-pixel particle (10-pixel equivalent diameter) represents a relative error of ~1%, manually traced images were processed at a resolution that ensured all plagioclase crystals with L > 0.15 mm were >80 pixels in area. Similarly, QEMSCAN images were collected at a resolution of 10 μm to ensure that plagioclase crystals with L > 0.15 mm and realistic length/width aspect ratios of three cover ≥75 pixels. Although QEMSCAN images could have been collected at higher resolutions, much longer acquisition times would have been required; imaging at a resolution of 2 μm would have taken ~25 times longer than imaging at a resolution of 10 μm.

The classic presentation of CSD data on semi-logarithmic plots of ln[n(L)] vs. L requires stereological conversion of 2D crystal size data. In order for such conversions to be accurate, crystals must define a single population of known and constant morphology (Higgins 2000). Although it is possible to make reasonable assumptions about crystal morphology within a single population based on crystal length and width data (Morgan and Jerram 2006), such assumptions cannot be applied across the multiple plagioclase populations present in the Laki lava (Guilbaud et al. 2007; Neave et al. 2013).

To identify different macrocryst populations while avoiding the pitfalls of stereological conversion, CSDs determined from all images were evaluated first on semi-logarithmic plots of number area density (N_A) normalized by bin width (bw) vs. the square root of crystal area (A^0.5). By plotting crystal size as A^0.5 rather than L, the effect of crystal morphology on stereologically uncorrected CSDs can be reduced (Neave et al. 2014). Furthermore, the following two binning strategies were used to ensure that our interpretations were not affected by how we chose to present the data: linear binning with a spacing of 0.05 mm (Armienti 2008), and geometric binning where each successive bin was a factor of 10^0.1 larger than the last (Sahagian and Proussevitch 1998).

Having identified coherent plagioclase populations using stereologically uncorrected CSDs, we then calculated classic CSDs (ln[n(L)] vs. L) from traced thin section images. Best-fitting crystal shapes were first estimated from corrected plagioclase lengths and widths using CSDslice (Morgan and Jerram 2006). These best-fitting shapes were then used with corrected plagioclase lengths to calculate geometrically binned CSDs using CSDCorrections (Higgins 2000). Owing to the small number of macrocrysts present (<75 with L > 0.5 mm after conversion), stereological conversions were only applied to statistically robust microcryst populations.

Total crystal contents and phase proportions determined by manual point-counting and QEMSCAN imaging are shown in Figure 4. Total crystal contents estimated by manual point-counting vary between 8.4 and 28.4% by area, and correlate loosely with surface transport distance and position in the eruption chronology. In contrast, total crystal contents estimated from QEMSCAN imaging vary between 32.9 and 64.4% by area. Phase proportion estimates also differ significantly between manual and automated techniques, with QEMSCAN data sets returning more plagioclase than manually obtained data sets (Fig. 4). Reasons for these discrepancies are discussed below. Total crystal contents and phase proportions estimated from QEMSCAN analyses performed at different resolutions (10–200 μm) are indistinguishable (Fig. 4).

Independent estimates of plagioclase contents from traced thin section images are compared with estimates from manual point-counting and QEMSCAN analyses in Figure 5. Plagioclase proportions determined by manual point-counting and image tracing generally agree; with the exception of plagioclase proportions from macrocryst-poor LAK09, estimates from image tracing are within the 2σ uncertainty of estimates from manual point-counting. However, plagioclase proportions determined by QEMSCAN imaging exceed those from manual methods by factors of 2.2–12.0. For example, while manual point-counting and image tracing return plagioclase content estimates of 6.5 and 6.4% for LAK04, respectively, the Q10 image contains 22.8% plagioclase by area. Granulating this Q10 image and discarding plagioclase particles with L < 0.15 mm, which would exclude all groundmass grains if segmentation were perfectly efficient, leads to a reduction in estimated plagioclase proportions to 18.0% by area (Q10G; Fig. 5). Manual rectification of granulated images to break apart glomerocrysts before discarding particles with L < 0.15 mm reduces estimated plagioclase proportions to 14.1% by area (Q10R; Fig. 5). However, even after manual rectification, plagioclase contents estimated by automated methods still exceed estimates from manual techniques in all samples by factors of 1.6–4.3.

Processing QEMSCAN images with SIP files tuned to distinguish between micro- and macrocryst, and groundmass plagioclase compositions (>An₆₅ and <An₆₅, respectively) leads to significant improvements in plagioclase phase proportion estimates (Q10S; Figs. 5 and 6). After tuning SIP files, estimated plagioclase contents are 8.1–14.0% by area, factors of 0.8–2.7 different from manual estimates. However, several falsely identified high-anorthite pixels still occur within sample ground-masses, which is consistent with the tendency to overestimate plagioclase contents (Fig. 6).

Plagioclase aspect ratios (length/width) determined from traced images are plotted against crystal size (A^0.5) in Figure 7. Aspect ratios are not shown for QEMSCAN data sets because segmentation was insufficiently reliable. Plagioclase morphology varies systematically as a function of grain size: macrocrysts (A^0.5 > 0.2 mm, which is equivalent to L > 0.5 mm) are equant and have mean aspect ratios of 2–3, whereas microcrysts (A^0.5 < 0.2 mm, which is equivalent to L < 0.5 mm) are elongate and usually have aspect ratios >4. Given that the shape of plagioclase grains reflects their crystallization histories (Lofgren 1974; Higgins 1996; Holness 2014), these size-morphology relationships confirm the presence of multiple plagioclase populations (cf. Neave et al. 2013).

QEMSCAN imaging combines textural and compositional information within single data sets, a notable advantage over other textural techniques. To characterize plagioclase size-composition relationships, we plotted the mean anorthite content of plagioclase particles as a function of their size (Fig. 8). Anorthite contents were estimated by calibrating granulated and rectified QEMSCAN images against the compositions used to generate SIP file entries. However, given the small number of plagioclase species that can be distinguished using the low count spectra necessary for the timely production of QEMSCAN images, these calibrations are only semi-quantitative. Robust size-composition relationships can nonetheless be identified in all samples: large particles are dominated by high-anorthite contents whereas small particles have variable but, on average, lower anorthite contents.

Porphyritic samples (LAK18, LAK09, and LAK04) have similar size-composition systematics (Figs. 8a–8c). The largest macrocrysts (A^0.5 > 0.4 mm) have highly anorthitic mean compositions (An₈₄–An₉₂) that lie within the range of published high-anorthite core compositions for Laki (Guilbaud et al. 2007; Neave et al. 2013). The large scatter in apparent microcryst compositions (A^0.5 < 0.2 mm) notwithstanding, the mean anorthite contents of these grains are slightly higher than those reported previously for microcrysts and macrocryst rims (An₆₇–An₇₀ vs. <An₆₅; Guilbaud et al. 2007; Neave et al. 2013). This offset probably relates to the frequent false identification of primitive bytownite pixels (~An₈₀) in zones of labradoritic composition (~An₆₀; Fig. 6a). Seriate LAK27 shows a similar trend in plagioclase size-composition space to the other samples that is offset to systematically lower anorthite contents; all plagioclase particles contain a greater proportion of labradorite pixels in LAK27 (Fig. 8d), reflecting the extensive groundmass crystallization experienced by this sample.

Plagioclase size distributions calculated from traced images are plotted in Figure 9 using diamond symbols. Most plagioclase crystals lie within the size range A^0.5 = 0.1–1 mm, though a few larger grains are present in LAK27. All CSDs are kinked, with microcryst and macrocryst populations defining distinct arrays. The kink separating these populations occurs at A^0.5 ~ 0.3 mm in most samples but is shifted to A^0.5 ~ 0.4 mm in LAK27 because of this sample's higher crystallinity (Fig. 5; Cashman and Marsh 1988; Higgins 1996). Note, however, that mean microcryst sizes are significantly smaller than the sizes indicated by kink points; kink positions are not equivalent to the size thresholds used to distinguish between plagioclase populations (L = 0.5 mm; A^0.5 = 0.2 mm). Grains smaller than A^0.5 ~0.1 mm are below the L = 0.15 mm threshold applied during tracing, resulting in an artificial drop in population density at the smallest grain sizes. CSDs calculated from linearly and geometrically binned data show similar trends for all samples; fits to plagioclase populations are independent of binning style when data are sufficiently dense. Kinks between microcryst and macrocryst populations also occur at similar A^0.5 values in the differently binned data sets. For geometrically binned data, intercepts of fits to microcryst populations (i.e., values of N_A/bw at A^0.5 = 0 mm), which are proxies for nucleation densities (n⁰), lie between 148 and 403 mm^–3. Intercepts of fits to macrocryst populations lie between 2.7 and 12.2 mm^–3.

Plagioclase size distributions calculated from unrectified and rectified Q10 images are plotted in Figure 9 using square and circular symbols, respectively. QEMSCAN-derived CSDs show the same primary features as CSDs calculated from traced thin section scans: populations of small and large plagioclase crystals separated by breaks in CSD slopes at A^0.5 ~0.3 mm. However, there are two important differences between CSDs obtained from tracing and from QEMSCAN images. First, plagioclase number densities estimated from Q10 images are substantially higher than those from traced images. For example, intercepts of fits to microcrysts in both unrectified and rectified Q10 CSDs lie between 1097 and 2981 mm^–3 in comparison with values of between 148 and 403 mm^–3 from CSDs from traced images. And, second, CSDs from QEMSCAN images are more smoothly concave than traced CSDs. Reasons for these differences are discussed below.

Plagioclase morphology varies as function of plagioclase size in the Laki lava (Fig. 7). Therefore, a single stereological correction cannot not be applied across the full range of plagioclase sizes present in each sample: currently available stereological correction schemes assume a uniform shape across all crystal sizes (Higgins 2000). Classic CSDs for estimating timescales of magmatic processes were thus calculated using plagioclase morphologies estimated from microcrysts (L < 0.5 mm). Separate conversions for macrocryst populations could not be performed because our samples do not contain sufficient numbers of macro-crysts to reconstruct 3D grain shapes robustly (>75; Morgan and Jerram 2006). Furthermore, macrocrysts are more primitive than microcrysts (Fig. 8) and record events that occurred before the final assembly and eruption of the Laki magma; chemical and structural complexity within large plagioclase macrocrysts from the Laki lava has been interpreted in terms of magma mixing and crystal mush entrainment (Guilbaud et al. 2007; Passmore et al. 2012; Neave et al. 2013). It is thus unclear what geological meaning could be extracted from the apparent crystallization timescales of large crystals that may have been resident in the plumbing system for thousands of years (e.g., Cooper et al. 2016).

Stereologically corrected CSDs are shown in Figure 10. Linear fits through microcryst populations in the porphyritic samples (LAK18, LAK09, and LAK04) have similar intercepts (n⁰ = 394–460 mm^–4), whereas the intercept for seriate LAK27 is slightly higher (n⁰ = 667 mm^–4). Gradients (–1/Gτ) of fits through data from episodes I and III (LAK18 and LAK09) are somewhat more negative than the gradients of fits through data from episodes VII and VIII (LAK04 and LAK27): –1/Gτ = –12.2 to –11.1 mm^–1 vs. –9.5 to –8.4 mm^–1, respectively.

Total crystal contents and phase proportions determined by manual and automated point-counting techniques differ significantly (Figs. 4 and 5). We suggest two main reasons for this: first, microcrysts were not discriminated efficiently from groundmass grains by automated point-counting; and second, automated processing methods are poor at segmenting glomerocrysts.

A minimum size threshold of L = 0.15 mm was imposed when distinguishing microcrysts from groundmass grains during manual point-counting (Passmore et al. 2012). However, ground-mass grains were not distinguished from microcrysts during the initial processing of QEMSCAN images; all pixels identified as plagioclase, clinopyroxene, or olivine contributed toward estimates of total crystal content. While points identified manually as groundmass encompassed glass, mesostasis, and fine-grained groundmass crystals, only glass and mesostasis were categorized as groundmass during QEMSCAN imaging.

Separating micro- and macrocryst pixels from groundmass pixels based on their composition was the most successful method for bringing manual and automated data sets into alignment: tuning SIP files resulted in plagioclase contents 0.8–2.7 times the magnitude of contents estimated using manual techniques. Numerous misclassified pixels (e.g., low-anorthite pixels in macrocryst cores) nonetheless remained, indicating that higher count spectra (>5000 counts) would be required to discriminate between compositionally different plagioclase populations robustly. Alternatively, QEMSCAN phase maps could be converted into grain maps by integrating crystal orientation information from electron backscattered diffraction (EBSD) analyses (e.g., Prior et al. 1999; Cordier et al. 2014).

Population densities determined from traced and QEMSCAN images are compared in Figure 11. Plagioclase population densities calculated from QEMSCAN images are almost always higher than those calculated from traced images regardless of the binning style or degree of rectification (Figs. 11a and 11b). The difference in population density between manual and automated methods also varies as a function of plagioclase size (Figs. 11c and 11d): QEMSCAN population densities [ln(N_A/bw)] are up to three log units higher at small grain sizes (A^0.5 ~ 0.1 mm), but only half a log unit higher at larger grain sizes (A^0.5 > 0.2 mm).

Although CSDs from manually traced images are susceptible to some uncertainties (e.g., tracing precision and imaging resolution), the linearity of microcryst CSDs indicates that no size-dependent sampling biases are present over the A^0.5 = 0.1–0.3 mm interval. In contrast, the inefficient segmentation of small plagioclase particles in QEMSCAN images is especially notable at the smallest grain sizes in which an abundance of groundmass plagioclase agglomerations biases calculations and results in concave-up CSDs. Therefore, CSDs from QEMSCAN images are most comparable with their manual counterparts at larger grain sizes (A^0.5 > 0.2 mm) where both automated segmentation and manual rectification processes are more reliable.

CSDs calculated from QEMSCAN images of porphyritic rocks are thus prone to greater uncertainties than those calculated from traced images because it is challenging to produce grain maps automatically. Nonetheless, even CSDs calculated from unrectified QEMSCAN images reproduce the samples' two most important textural features: the division of plagioclase crystals into microcryst and macrocryst populations, and the higher number density of macrocrysts in samples from episodes VII and VIII than in samples from episodes I and III. Although CSDs from labor-intensive, manually collected data sets are still required for accurate quantification, CSDs from automatically collected data sets can nevertheless provide a rapid overview of textural features with minimal user input.

Although the textural properties of volcanic rocks can be used to investigate deep magmatic processes (Higgins and Roberge 2007; Armienti et al. 2013; Fornaciai et al. 2015), it is also important to assess the effects of surface transport on the evolution of rock textures. This is especially important in the case of the Laki eruption because thermally efficient transport over long distances resulted in continued evolution of the lava during emplacement (Guilbaud et al. 2007). For example, combining the minimum rate of lava surge advance (2 km per day) with a maximum lava transport distance (~40 km) indicates that some lava batches may have been transported within channels and lava tubes for at least 20 days before final emplacement (Thordarson and Self 1993). Indeed, surface transport timescales on the order of days are supported by timescales of diffusive H₂O loss from olivine-hosted melt inclusions collected from lava selvages 5–30 km downstream from the eruption site (4.0 ± 3.4(1σ) days; Hartley et al. 2015).

Samples were collected at different distances from their source vents (Fig. 1), implying that they experienced different post-eruptive thermal histories. Although plagioclase contents roughly correlate with transport distances (LAK09, which was collected at <1 km from its source vent, has the lowest crystal content, whereas LAK27, which was collected at ~40 km, has the highest), most microcrysts are still too large to have grown during emplacement and must reflect earlier magmatic processes. However, a mixture of pre- and post-eruptive processes could be recorded in the size distributions of the plagioclase microcrysts (L = 0.15–0.50 mm).

For stereologically corrected CSDs (Fig. 10), gradients of fits through plagioclase microcryst populations (–1/Gτ) become shallower with increasing distance from the samples' source vents: –12.2 mm^–1 in LAK09 (<1 km); –11.1 mm^–1 in LAK18 (~5 km); –9.5 mm^–1 in LAK04 (~15 km); and –8.4 mm^–1 in LAK27 (~40 km). Assuming that all samples contained similar microcryst populations at the time of eruption, as implied by the compositional homogeneity of microcrysts ejected throughout the course of the eruption (Guilbaud et al. 2007; Neave et al. 2013), these textural differences are best accounted for by varying degrees of CSD ripening during surface transport (Cashman and Marsh 1988; Higgins and Roberge 2003): LAK04 and LAK27 were transported further prior to final emplacement than LAK18 and LAK09, and thus experienced a greater degree of post-eruptive textural modification.

The relatively albitic composition of plagioclase microcrysts (~An₇₀; L < 0.5 mm) reflects their growth from liquids closely associated with erupted tephra glasses (Guilbaud et al. 2007; Neave et al. 2013). In contrast, macrocrysts (L > 0.5 mm) are dominantly composed of high-anorthite cores (An₈₄–An₉₂; Fig. 8) that are far from being in equilibrium with tephra glasses. Indeed, these high-anorthite compositions cannot be related to the erupted glasses by evolution along a simple single liquid line of descent (Neave et al. 2013): no melts with sufficiently high Ca/Na values to stabilize high-anorthite plagioclase can be generated by adding observed phase compositions back into the magmatic tephra compositions. By combining these findings with records of trace element heterogeneity in olivine-hosted melt inclusions associated with high-anorthite plagioclase cores, Neave et al. (2013) suggested that primitive macrocrysts were sourced from disaggregated crystal mushes formed from a distribution of mantle melts different from that which formed the erupted liquid. The CSDs and size-composition information (Figs. 8–10) presented here provide strong independent evidence for the pre-eruptive entrainment of disequilibrium macrocrysts into the Laki carrier liquid.

Although episodes I–V were the most vigorous and productive of the Laki eruption (Thordarson and Self 1993), lavas from these episodes contain the fewest crystals (<9% by area; Fig. 14b in Passmore et al. 2012). In contrast, lavas from fissures VII–VIII carry considerable crystal loads (mean 19% by area; Fig. 14b in Passmore et al. 2012). Fits through macrocryst populations in CSDs from LAK04 and LAK27 (fissures VII and VIII) have shallower gradients and higher intercept values than similar fits through macrocryst populations in LAK18 and LAK09 (fissures I and III) (Figs. 9 and 10), implying that macrocrysts are not only larger in the products of these later episodes, but that they are also more abundant. Thus, given that only modest differences between microcryst populations from different samples can be attributed to post-eruptive processes, the total crystal content of lava samples reflects primarily the abundance of macrocrysts and hence the amount of mush entrainment at depth.

When combined with the full point-counting data set of Passmore et al. (2012), our CSDs suggest that the efficiency of mush disaggregation and consequent macrocryst entrainment increased significantly between episodes I–III and episodes VII–VIII. A lower macrocryst content in episodes IX and X also suggests that macrocryst entrainment efficiency may have decreased toward the end of the eruption (Passmore et al. 2012). While variability in disaggregation efficiency may represent lateral variations in mush and cumulate petrology at depth, dynamic processes may have also been important. However, evaluating the controls on mush disaggregation efficiency is beyond the scope of this study, and may be approached better using numerical approaches (e.g., Bergantz et al. 2015; Schleicher et al. 2016).

Plagioclase microcrysts record the final phase of crystallization before eruption that, according to crystal zoning patterns, must have occurred immediately after the entrainment of large, primitive macrocrysts. Stereologically corrected CSDs of plagioclase microcrysts thus record information about how timescales between disaggregation and ejection at the surface evolved as the eruption proceeded.

Once uncertainties in stereological conversions and the identification of coherent crystal populations have been minimized, crystal growth rates represent the largest source of error when estimating crystallization times from CSDs (Fornaciai et al. 2015). Plagioclase growth rates estimated from recent isobaric and isothermal crystallization experiments on mafic compositions are within the range G = 10^–8–10^–5 mm/s (Conte et al. 2006; Orlando et al. 2008; Agostini et al. 2013; Shea and Hammer 2013), with experiments carried out on an anhydrous trachybasalt from Etna at low degrees of under-cooling probably representing the most relevant growth rates for the Laki system (G ~ 0.5 × 10^–7 mm/s; Orlando et al. 2008). These growth rates are also broadly consistent with experiments carried out on a hydrous high-K basalt from Stromboli at low degrees of undercooling (G ~10^–7–10^–6 mm/s; Agostini et al. 2013). Thus, assuming 0.5–5 ×10^–7 mm/s as a range of feasible growth rates, we estimate that plagioclase microcrysts in LAK18 and LAK09 crystallized over 2–20 days. Equivalent crystals in LAK04 and LAK27 record timescales of 2.4–24 and 2.8–28 days, respectively. Differences in timescales between proximal (LAK18 and LAK09) and distal (LAK04 and LAK27) samples of up 0.4–8 days are however consistent with the degree textural of ripening expected to occur in the distal samples as a result of surface transport (see above). Thus, once the effects of post-eruptive modification have been accounted for, plagioclase microcrysts in all samples record coherent crystallization timescales of ~2–20 days, regardless of their position in the eruption's chronology.

Plagioclase microcryst crystallization timescales are equivalent to the timescales of rim formation on primitive macrocrysts and thus constrain their entrainment timescales (Neave et al. 2013). It is therefore encouraging that our 2–20 days residence timescale for primitive macrocrysts in the Laki carrier liquid brackets the 6–10 days timescale estimated by modeling the diffusive re-equilibration of olivine macrocryst rims from episodes I, III, and V (Hartley et al. 2016). Further validation of our 2–20 days disaggregation timescale estimate, and hence our choice of plagioclase growth rates, is provided by the elevated water content of primitive olivine-hosted inclusions in magmatic tephra samples: the extent of diffusive over-hydration observed in these inclusions requires that entrained macrocrysts spent a minimum of 2.5–19.1 days in the evolved Laki carrier liquid prior to eruption (Hartley et al. 2015). Moreover, the high mean aspect ratios of plagioclase microcrysts (length/width > 4) are indicative of crystallization timescales in the order of days to tens of days, whereas the low mean aspect ratios (length/width = 2–3) of large plagioclase macrocrysts are consistent with much longer crystallization timescales of years to hundreds of years, assuming that crystallization was continuous (Holness 2014). As noted by Hartley et al. (2016), a 2–20 days timescale is not only much shorter than the total duration of the eruption (245 days), but is also comparable in length with the intervals between eruptive episodes (1–28 days; Thordarson and Self 1993).

When compared with manual point-counting approaches, automated approaches can greatly overestimate the erupted crystal content of porphyritic samples unless phase identification algorithms are specifically tuned to take account of compositional differences between macrocrysts, microcrysts, and groundmass grains. Differences between CSDs generated by manual and automated techniques reflect the inefficient segmentation of glomerocrysts when processing QEMSCAN images, and thus highlight the limited utility of this technique for studying syneruptive processes recorded by the smallest of crystals. However, CSDs derived from automatically generated images do recapitulate the main features of CSDs from manually traced images, meaning that magma reservoir processes recorded by larger macrocrysts are suitable for investigation with automated mineralogical methods. Many textural properties of volcanic rocks can thus be estimated with a fraction of the user input required for traditional methods. For example, key samples suitable for high-resolution but labor-intensive manual analysis can be selected from larger samples suites by performing automated analyses beforehand: while generating CSDs from thin section images can take tens of hours, generating CSDs from QEMSCAN data sets takes a matter of minutes (though collecting high-resolution QEMSCAN images still requires several hours of instrument time).

Large (L > 0.5 mm), low aspect ratio (length/width = 2–3) macrocrysts from the Laki lava flow contain uniformly primitive cores (An₈₄–An₉₂), whereas small (L < 0.5 mm), high aspect ratio (length/width >4) microcrysts are always more evolved (<An₈₀) and approach compositions in equilibrium with the erupted melt (~An₆₅). Large plagioclase macrocryst cores are too anorthitic to be related to the erupted melt by simple fractional crystallization and were probably sourced from disaggregating crystal mushes. Variations in macrocryst contents and CSDs between samples demonstrate that macrocryst entrainment efficiency varied during the course of the eruption: samples from later eruptive episodes (VII and VIII) carry a greater mush-derived component than samples from early episodes (I and III).

Compositional zoning patterns indicate that microcrysts grew immediately after the entrainment of macrocrysts and can therefore be used to constrain timescales between mush disaggregation and eruption. In turn, these timescales can be used to test whether mush disaggregation occurred in a single event before the eruption started (Scenario 1 in Fig. 12), or throughout the course of the eruption (Scenario 2 in Fig. 12). Using geologically plausible plagioclase growth rates, we estimate 2–20 days timescales for the simultaneous growth of microcrysts and entrainment of macro-crysts that are in good agreement with estimates from independent geochronometers. These timescales are also comparable with inter-episode repose times, implying that primitive macrocrysts erupted in later episodes were locked within mushes when the eruption started. We therefore conclude that mush disaggregation occurred throughout the course of the eruption (Scenario 2 in Fig. 12), with magmas from each episode being assembled on the order of 10 days before being ejected from their respective source vents. Importantly, an approximately 10-day time frame corresponds with historical records of seismicity before the onset of many eruptive episodes (Fig. 12; Thordarson and Self 1993, and references therein). Pre-eruptive seismicity may have conceivably been generated by the same magma movements that resulted in mush disaggregation; petrological and geophysical expressions of magmatism appear to be related at Laki.

Our textural observations imply that the near-homogenous composition of the Laki magma cannot have formed during a single mixing event (Scenario 1 in Fig. 12). With the sole exception of mush entrainment efficiency, magma assembly and evolution processes must have thus remained remarkably consistent throughout the eruption because the products of different eruptive episodes have very similar compositions despite their separation in space and time (Passmore et al. 2012). That is, successive batches of erupted magma must have crystallized from similar parental magma distributions under similar P-T conditions, suggesting that there were no substantial changes in reservoir architecture over an eight-month period.

Our findings also suggest that the magma feeding the 15.1 km³ eruption was mobilized in a punctuated manner. Specifically, the short entrainment and mush disaggregation timescales we calculate imply that magma batches from each episode were only mobilized a matter of days before their eruption. Indeed, close temporal relationships between eruptions and ground deformation events during the 1977–1984 CE Krafla Fires in north Iceland demonstrate that repeated phases of magma movement and, conceivably, mush entrainment are unlikely to be unique to Laki (Björnsson 1985). Moreover, compositional heterogeneity in the products of the 1730–1736 CE Timanfaya eruption on Lanzarote and 871 and 1477 CE Veiðivötn eruptions in Iceland confirm that fissure eruptions are often fed from plumbing systems in which communication between different magma batches is incomplete or completely absent (Carracedo et al. 1992; Sigmarsson et al. 1998; Zellmer et al. 2008).

One key implication of the progressive mobilization we infer here is that the volume of eruptible magma in the Laki plumbing system at any given time was probably much closer the 1–2 km³ erupted per episode than the final erupted volume of 15.1 km³. Therefore, each episode probably involved a volume of mobilized magma comparable to that which fed the 1.5 km³ 2014–2015 Bárðarbunga-Holuhraun eruption (Guðmundsson et al. 2016). Given that Laki's eruptive episodes were similar in size to numerous documented eruptions from Iceland and elsewhere, the eruption's exceptionalism in the recent geological record may thus hinge more on its tremendous vigor than on the ultimate volume of its products: during some episodes, the same volume of lava was emplaced within in a few weeks as was emplaced during the whole six-month duration of the 2014–2015 Bárðarbunga-Holuhraun eruption.

This work was supported by NERC grants NER/S/A/2004/12727, NE/1528277/1, and NE/I012508/1. D.A.N. acknowledges support from the Alexander von Humboldt Foundation. We thank Margaret Hartley for her comments on an earlier draft of this manuscript and discussions on the Laki eruption, Simone Tommasini and one anonymous reviewer for their detailed and constructive comments, and Chiara Maria Petrone for her helpful guidance and efficient editorial handling.