Forecasting the timing and size of volcanic eruptions requires a proper interpretation of multiparametric monitoring signals. Studies of the erupted rocks can provide critical information on the processes and volcano plumbing system that is needed to decode the monitoring signals. Here we present the results of a petrological study of plagioclase phenocrysts using a new statistical approach that allows us to estimate the amount of intruded magma before eruption. Our crystal population analysis of the 2006 and 2010 CE Merapi volcano (Indonesia) eruptions shows that ∼60 ± 20 vol% of the 2010 magma was left over from the 2006 magma, and thus ∼40 ± 20 vol% was newly intruded magma. Using the published values of the 2010 erupted magma volume, this corresponds to >8 to 20 (±7) × 106 m3 of new magma. This is a minimum estimate and is similar to the inferred pre-eruptive deformation volume (18 ×106 m3), although given the uncertainties, several million cubic meters of magma intruded in 2010 could still be in the Merapi plumbing system. Our approach could be used at other volcanoes to quantify the volume of intruded magma and thus help in better understanding the unrest signals that anticipate eruptions.

Many volcanoes erupt magmas of the same composition for decades or centuries, e.g., Mayon in the Philippines (Newhall, 1979), Merapi in Indonesia (Costa et al., 2013), Arenal in Costa Rica (Reagan et al., 1987), or Soufrière Hills in Montserrat (Murphy et al., 2000). What drives these frequent eruptions of the same magma? How much of the magma produced by each eruption is new replenishment, and how much was already stored in the system? Being able to answer these questions is important for understanding how volcanos work, and thus for being able to properly interpret monitoring data and unrest signals.

Geophysical studies can provide quantitative constrains on the amount of intruded magma stored before eruption via inverse modeling of the deformation field (e.g., Segall, 2010) and from the cumulative seismic moment (e.g., White and McCausland, 2016). Changes in the gas flux and mass balance can in principle also provide information on the amount of resident and replenishment magma (e.g., Allard et al., 1994). However, many active volcanoes that erupt frequently do not show significant deformation or seismicity before eruption (e.g., Biggs et al., 2014). Moreover, magma volume-change estimates before eruption can be nonunique because of the tradeoff between various parameters and the properties of the host rock (Segall, 2010).

An alternative source of information to constrain the amount of intruded magma comes from detailed petrological and geochemical study of erupted rocks (e.g., Druitt et al., 2012; Humphreys et al., 2013). The identification of the various magma components is relatively straightforward if the resident and intruding magma are very different (e.g., quenched mafic inclusions in silicic host), including having very different isotopic signatures (e.g., Alloway et al., 2015). But many volcanoes erupt magmas with very similar petrological and geochemical features, which makes quantitative estimation of resident and intruded magma a very difficult task. Here we show that by analyzing a large number of plagioclase phenocrysts and applying a newly developed tool based on statistical analysis, it is possible to identify the volumes of resident and intruding magmas. Such information can be then combined with seismic or geodetic estimates to better interpret monitoring signals and thus better anticipate eruptions.

Background on the Merapi 2006 and 2010 CE Eruptions

Merapi volcano (Java, Indonesia) is one of the most active and hazardous volcanoes in the world. It is known for its small to moderate-size eruptions that alternate between effusive and explosive every 4–6 yr (Surono et al., 2012). The 2006 eruption was mainly effusive and relatively small (∼5 × 106 m3 of dense rock equivalent [DRE]; Ratdomopurbo et al., 2013), in contrast to the much larger (5–20 × 106 m3 DRE) and explosive eruption of 2010 (Surono et al., 2012). Such a difference in volume and eruption style is, however, not reflected in a major change in bulk-rock compositions, including isotopes (Handley et al., 2014), most minerals (Costa et al., 2013), or pre-eruptive conditions (e.g., Costa et al., 2013; Erdmann et al., 2016). Petrological evidence for new magma replenishment in 2010 includes dissolution textures in plagioclase and clinopyroxene, the presence of Al-rich amphibole (among other observations) (Costa et al., 2013; Erdmann et al., 2014), and the more primitive magma of the last phase of the eruption (Komorowski et al., 2013; Erdmann et al., 2016).

Samples and Analysis of Plagioclase Composition

We used two samples from the 2006 eruption: 1506L2, which we refer to as 2006A; and 1406L2, which we label 2006B (Costa et al., 2013). Both samples were collected in Kali Gendol and erupted on two different days in June 2006. We also used two samples from the 2010 eruption: 1010523, which we label 2010A; and 101121-2B, which we label 2010B (Costa et al., 2013). Sample 2010B is likely from the dome that formed during 1–4 November (eruptive stage 3 of Komorowski et al. [2013]). Sample 2010A is a bread-crusted bomb in the pyroclastic flow from the main explosive phase of 5 November 2010 (eruptive stages 4–6 of Komorowski et al. [2013]). These 2010 samples were previously studied by Costa et al. (2013), and we believe that they are representative of the magma that formed the largest volume of the eruption, although not of the last eruptive phase that contains the more primitive and crystal-rich magma (stage 7; Komorowski et al., 2013; Erdmann et al., 2016).

We obtained back-scattered electron (BSE) images of four full petrographic thin sections using a scanning electron microscope (Nanyang Technological University, Singapore; NTU) (see the GSA Data Repository1 for more detail). BSE images of plagioclase were calibrated for anorthite content (An = 100 × Ca/[Ca + Na] in mol) with multiple individual analyses by electron microprobe (at NTU). We extracted plagioclase phenocryst sections (crystals >0.1 mm) from BSE images using CEmin software (Zeng et al., 2018). In total, we obtained 243 and 236 plagioclase phenocryst compositional data for each sample of the 2006 eruption respectively, and 273 and 283 phenocrysts for each sample of the 2010 eruption. These numbers are high enough to obtain representative sampling of the crystal populations (Cheng et al., 2017).

Identification and Classification of Plagioclase Populations

To identify the various plagioclase populations, we used a recently developed approach by Cheng et al. (2017). The most complex part of this approach is to filter out the effect of random cuts of three-dimensional (3-D) plagioclase crystals into two-dimensional (2-D) sections from the real variety caused by the presence of multiple populations. First, we acquired a large number of An histograms of plagioclase sections from a given sample. Then we computed the aggregate An distribution for the whole sample and calculated the similarity between the distributions of the overall sample and of each plagioclase section (expressed as the normalized difference between two An histograms). We identified the sections that had the highest similarities and then applied thresholds to find the main peaks of similarity that are representative of the sample. Such peaks allow the identification of the so-called reference (REF) and ideal (ID) sections, which correspond to the real 3-D crystal populations (Cheng et al., 2017). Additional details on the method are found in the Data Repository and in Cheng et al. (2017).

Volume Calculations

We used the proportions of the same plagioclase population (e.g., Pop1) from two consecutive eruptions (e1 and e2; i.e., 2006 and 2010) and the volume of the last eruption (Ve2; i.e., 2010) to calculate the minimum volumes of mingled magmas (resident and new intrusion). The basic relation is that Vr / (Vr + Vin) = [%Pop1 (in e2)] / [%Pop1 (in e1)], where Vr is the volume of resident magma and Vin is the volume of intruded magma. These volumes are unknown, but the percent of each plagioclase population (e.g., %Pop1 in e1) that is shared between the two magmas is measured from the erupted rocks. For simplicity, we assume that Pop1 belongs exclusively e1 (but this is not a necessary condition of our model), which gives the minimum value of Vin. The volume of erupted magma in e2 that is recycled from the resident magma is Ve2 × [%Pop1 (in e2) / %Pop1 (in e1)], and that of intruded magma is Vin = Ve2Vr (see the Data Repository).

Sources of Uncertainty

Uncertainties can arise from the petrological variability of the studied samples. From image analysis, we obtained the following plagioclase phenocryst volume proportions: 22% for sample 2006A and 29% for 2006B, and 22% for sample 2010A and 20% for 2010B. Thus, there is a variability of ± 4 vol% for the 2006 eruption, and ± 1 vol% for the 2010 eruption. The uncertainty on the histograms depends on the total counts and pixel size of BSE images, and on the calibration of the grayscale on the An composition maps, which together give errors of ∼1% (Cheng et al., 2017). The uncertainty in the identification of the various populations depends on the sample itself. Cheng et al. (2017) found that if there is only a single population, 10% or less of these crystals could be another population (i.e., identification is at least 90% correct), depending on the sample. When there are more populations, we give the standard deviation of the best-fit solution, which is typically 10%–20% of the ratios of the populations. We obtained an uncertainty of 11% for the ratio of two populations of sample 2010A and 15% for ratio of two populations of sample 2010B. There is also variability in the amount of plagioclase and in different samples.

The plagioclase phenocryst sections of the 2006 and 2010 samples show a wide range of zoning patterns, with numerous resorption surfaces, shifts in An content, and geometries (Fig. 1; see the Data Repository). A large part of the compositional and textural variety of the 2-D sections is apparent, and is mainly caused by the random sampling of 3-D plagioclase crystals with complex geometries and zoning patterns as seen in 2-D sections (Cheng et al., 2017). The histograms of An of all aggregate plagioclase sections for each of the two samples from the 2006 eruption (2006A and 2006B) are 92% similar, and both have close to a normal distribution with a mode at approximately An53 (An = 53; Fig. 2A). The An histograms of the two samples from the 2010 eruption (2010A and 2010B) also have a high similarity (89%), but with a narrower distribution than the 2006 samples and a mode at above An50 (Fig. 2). The similarity between the overall An histograms of the 2006 and 2010 deposits is very high, ranging from 80% to 88% (85% ± 3.5%; Table 1).

The statistical similarity analysis between all of the plagioclase sections of the two samples of the 2006 deposits shows that 91% of the An distribution of sample 2006A and 96% of the An distribution of sample 2006B can be explained by the same plagioclase population (the so-called 2006A-Pop1 and 2006B-Pop1) with a very high similarity (90%; Fig. 2; Table 1). The remaining 4%–9% of phenocrysts could belong to other low-abundance populations that we can’t identify with our approach.

The similarity analysis of the 2010 samples shows that they each have two plagioclase populations: for sample 2010A, 2010A-Pop1 makes up 71% ± 11% of the data, and 2010A-Pop2 makes up the remaining 29 vol%; for sample 2010B, 2010B-Pop1 makes up 70% ± 15% of the data, and 2010B-Pop2 makes up 30% (Fig. 2). Moreover, for both 2010 samples, we can reproduce the An distribution of 93%–94% of the crystals, so that ∼7% could belong to other unidentified low-abundance populations (2010-PopOther). The similarities between the Pop1 of the two samples and between the Pop2 of the two samples is high (83%–84%; see Table 1 for full similarity table). By definition, the similarity between the two different populations (Pop1 and Pop2) is much lower (66% and 58%; Table 1).

We quantified the relation between the 2006 and 2010 magmas by calculating the similarities of their plagioclase populations (Table 1). We found that Pop1 from 2006 and Pop1 from 2010 deposits are quite similar, with similarity values of 81% ± 2%, which we take as indication that they are basically the same crystal population. This suggests that Pop1 in 2010 is recycled from 2006 magma left in the system. Such crystal recycling is commonly reported from frequently active volcanoes at subduction zones (e.g., Jerram and Bryan, 2015). In contrast, the similarity between Pop1 of 2006 and Pop2 of 2010 is lower at 70% ± 3%, and we suggest that Pop2 from 2010 represents the new magma input. This implies that magma mingling was quite effective and may have been facilitated by heat addition and arrival of hot volatiles from primitive magma (e.g., Spera et al., 2016). This means that 70 vol% of the phenocrysts erupted in 2010 (those corresponding to Pop1) could have been recycled from leftover magma from 2006 in the plumbing system, and opens the door to estimating the amount of magma intruded in 2010.

The identification of the various plagioclase crystal populations in the 2006 and 2010 magmas allows us to propose a magma plumbing model for Merapi and the volume relations of the different components (Fig. 3). The estimated 2010 erupted magma (DRE) is 20–50 × 106 m3 (Surono et al., 2012). Analysis of the deformation patterns prior to the 2010 eruption (Aisyah et al., 2018; Kubanek et al., 2015) shows an increase in volume of 18 × 106 m3, which almost overlaps with the estimated range of erupted magma volume, but it could be up to a factor of 3 smaller. The difference between the erupted and deformation volumes is common to many other volcanoes and eruptions (e.g., Kilbride et al., 2016) and could be partly explained by the effect of magma compressibility.

Our petrological analysis based on plagioclase populations shows that 30% ± 13% of plagioclase is from newly intruded magma, the remaining being recycled leftover magma in the plumbing system from 2006. Based on this, it is possible to estimate the minimum volume of new magma (see the Methods section). The abundance of plagioclase in the 2010 samples is ∼21 (± 1) vol%, 70% (± 13%) of which is the resident magma, which corresponds to 15 ± 3 vol% of old plagioclase. Using the ratio of Pop1 in 2010 to that of Pop1 in 2006 as an indication of the mixing proportions between the resident and intruded magma, we find that the resident magma accounts for 60% ± 20% of the erupted magma volume in 2010 and 40% ± 20% was new magma.

Using these proportions and the erupted 2010 volume, we find that the newly intruded magma amounts to 8–20 (±7) × 106 m3 (Fig. 3). This volume estimate overlaps with that obtained from pre-eruptive deformation, although given the variations and errors in the various volume estimates, a few million cubic meters of the 2010 intrusion could still remain in the plumbing system. A recent tomographic study shows that seismic velocity changes below Merapi (Widiyantoro et al., 2018) could correspond to at least three magma reservoirs at different depths. The “2-D sizes” of the magma reservoirs are reasonably well constrained, but the deepest one appears to be several orders of magnitude larger than the erupted volume of 2010 (Widiyantoro et al., 2018). Thus, the deformation and petrological volume estimates we discuss are more likely representative of magmatic processes related to the two upper reservoirs: mingling probably happened in the shallower one (<4 km depth), and the new input may be from the intermediate reservoir at 10–20 km. Much more melt certainly remains below Merapi, from what we infer from our petrological analysis.

Our study shows that it is possible to combine geophysical and petrological data to better understand the magma budget and plumbing system of active volcanoes. This information provides better constraints for understanding and interpreting monitoring unrest signals that occur before eruptions. Our petrological estimates of magma volumes are especially useful for eruptions older than the monitoring records, and for volcanoes for which there are no deformation data or which don’t deform significantly (e.g., Biggs et al., 2014).

We are thankful to C. Widiwijayanti, and W. Li for many discussions about Merapi petrology and deformation. Reviews by L. Caricchi, S. Erdmann, three anonymous reviewers, and the editor (Chris Clark) were very helpful in clarifying various aspects of the manuscript. This research is supported by the National Research Foundation of Singapore and the Singapore Ministry of Education under the Research Centres of Excellence initiative, and a National Research Foundation of Singapore Investigatorship Award (NRF-NRFI2017-06) to F.C.

1GSA Data Repository item 2019398, details of the plagioclase similarity database and calculation procedures, is available online at http://www.geosociety.org/datarepository/2019/, or on request from [email protected].