Research ArticleGEOLOGY

Quantifying gas emissions from the “Millennium Eruption” of Paektu volcano, Democratic People’s Republic of Korea/China

See allHide authors and affiliations

Science Advances  30 Nov 2016:
Vol. 2, no. 11, e1600913
DOI: 10.1126/sciadv.1600913


Paektu volcano (Changbaishan) is a rhyolitic caldera that straddles the border between the Democratic People’s Republic of Korea and China. Its most recent large eruption was the Millennium Eruption (ME; 23 km3 dense rock equivalent) circa 946 CE, which resulted in the release of copious magmatic volatiles (H2O, CO2, sulfur, and halogens). Accurate quantification of volatile yield and composition is critical in assessing volcanogenic climate impacts but is challenging, particularly for events before the satellite era. We use a geochemical technique to quantify volatile composition and upper bounds to yields for the ME by examining trends in incompatible trace and volatile element concentrations in crystal-hosted melt inclusions. We estimate that the ME could have emitted as much as 45 Tg of S to the atmosphere. This is greater than the quantity of S released by the 1815 eruption of Tambora, which contributed to the “year without a summer.” Our maximum gas yield estimates place the ME among the strongest emitters of climate-forcing gases in the Common Era. However, ice cores from Greenland record only a relatively weak sulfate signal attributed to the ME. We suggest that other factors came into play in minimizing the glaciochemical signature. This paradoxical case in which high S emissions do not result in a strong glacial sulfate signal may present a way forward in building more generalized models for interpreting which volcanic eruptions have produced large climate impacts.

  • Volcanic gas emissions
  • volatiles
  • Paektu
  • Millennium Eruption
  • melt inclusions


Volcanic emissions have profound impacts on planetary atmospheres and drive climate change over a range of temporal and spatial scales (13). As the principal source of sulfur in the stratosphere (4), explosive eruptions that generate stratospheric clouds can result in long-lasting (in essence, several years) atmospheric effects via the injection of SO2, which oxidizes to form sulfate aerosol and promotes global cooling. Rhyolitic magmas that feed explosive volcanic eruptions, although characteristically S-poor, may contribute a substantial amount of S into the atmosphere sourced from a S-rich pre-eruptive vapor phase (513). Vapor saturation may be common in highly evolved silicic magmas that require protracted crystallization during crustal storage, resulting in the exsolution of vapor via second boiling (14, 15). Although the presence, amount, and composition of such pre-eruptive gas are not recorded directly in the rock record, it is possible, in some cases, to constrain it by examining geochemical trends in erupted samples (16, 17).

S yields and subsequent climate impacts from modern volcanic eruptions may be measured directly via satellite or in situ remote sensing. To evaluate the impacts of volcanoes on the atmosphere and climate throughout geologic time, we require a way to assess S yields of ancient or unmonitored eruptions. Many factors, including but not limited to gas yield and composition, dictate the degree of climate forcing imposed by a volcanic eruption. A large S yield may be required for, but is not necessarily indicative of, strong climate forcing. Other factors, such as latitude, the season of eruption, and the height of the eruption column, can determine the extent to which an S-rich volcanic eruption will perturb climate. Sulfate concentrations in polar ice cores and proper interpretation thereof may elucidate the efficiency of atmospheric transport of emitted volatiles and are thus thought to be the best records of changes in atmospheric chemistry through time. By combining independent analyses of (i) total gas yields from large explosive eruptions as recorded in volcanic rocks and (ii) perturbations in atmospheric chemistry, we can paint a more complete picture of the role of volcanism on climate.

Here, we demonstrate a method for constraining the total gas budget of large silicic eruptions by examining geochemical trends in glasses, crystals, and crystal-hosted melt inclusions (MIs), which represent various stages of the crystallization history of the magma. Petrologic and thermodynamic modeling reveals the evolution of volatiles in pre-eruptive melt and a coexisting vapor phase in the buildup to the comenditic Millennium Eruption (ME) of Paektu volcano, circa 946 CE. This technique considers gas generated before and during eruption plus contributions to the volatile budget from solid phases (for example, sulfide). Using this method, we calculate the total amount of gas generated during crustal magma storage, and so the values we report represent maximum possible volatile yields. We find evidence for the generation of a S-rich pre-eruptive vapor phase, which supplies most of the erupted S. If true gas yields were close to our maximum value, then our new gas yield estimates would place the ME among the largest emitters of climate-forcing gases in the Common Era (Fig. 1). This alone implies the potential for volcanogenic climate effects not unlike those seen after the 1815 eruption of Tambora (18), which was responsible for the “year without a summer” in 1816 and the deaths of more than 71,000 people. Conversely, however, relatively low-level sulfate deposits in Greenland ice cores recently linked to the ME (19) suggest only a modest S release to the atmosphere. Previous work citing low S yield estimates (2 Tg) appears consistent with this interpretation, but this value excludes the potential contribution from pre-eruptive fluid or solid-phase breakdown within ME magmas (20). Here, we demonstrate that a potential lack of strong climate forcing after the ME occurred in spite of the substantial S yield and suggest that other factors, such as latitude and seasonality, predominated in minimizing climatic effects.

Fig. 1 Map showing the location of Paektu volcano and associated tephra fallout and gas yields from large historic eruptions.

(A) Map showing the location of Paektu volcano (black triangle) and isopachs illustrating the extensive dispersal of tephra fallout from the ME [data from Machida et al. (61)]. Numbers within each oval indicate the thickness of tephra. Base map data from MapBox and OpenStreetMap. (B) Total S and halogen (F + Cl) yields from large historic eruptions, with eruption dates given in parentheses. Taupo, Katmai, Krakatau, and Tambora S yields from ice core (IC) (58). Pinatubo S yield from satellite remote sensing (RS) (72). Laki S and halogen yields estimated by the petrologic method over a period of 8 months (P) (73). Soufriere Hills and Mt. St. Helens (MSH) gas yields by the petrologic method (P) (74). The estimates for the Paektu ME are from this work (red bars) or by the petrologic method (blue bars) (20). Numbers above red bars indicate change in halogen and S yields estimated by this work compared to the estimates of Horn and Schmincke (20).

Volatile budgets of ancient silicic eruptions and the excess gas problem

Modern eruptions may be monitored and their gas yields measured directly via satellite or in situ remote sensing. For ancient or unmonitored eruptions, geochemists rely on clues in the rock record to determine the amount and composition of erupted gas. This is often achieved by comparing geochemical signatures in rocks that record pre- and post-eruptive volatile histories. Pre-eruptive magmatic volatile contents are commonly characterized using dissolved volatile concentrations in phenocryst-hosted MIs: small beads of liquid that become trapped within rigid crystal hosts and quench to glass upon eruption (21). Once enclosed within a crystal, the MI is largely isolated from changes in the magmatic system (in essence, magma differentiation and degassing) and so acts as a time capsule recording snapshots of melt chemistry, including dissolved volatiles, throughout the pre-eruptive evolution of the magma. Conversely, interstitial liquid between crystals [matrix glass (MG)] is privy to changes within the system and thus will degas most of its dissolved volatile content upon eruption as magma ascends and confining pressure decreases. The proportion of volatiles lost during eruption can thus be calculated as the difference between volatile concentrations in MG and those in the most evolved MI, which represent the magma just before eruption. If the volume and density of erupted material are known, then a mass of volatiles released during the eruption can be calculated. This is known as the “petrologic method” and is often applied to ancient or unmonitored eruptions for which no direct measurement of gas yield exists.

Advances in remote sensing (22) have provided a way to test the petrologic method by comparing calculated to measured gas yields. This has led to the discovery of what is known as the “excess S” or “excess gas” problem: S yields quantified by the petrologic method are often much lower than those measured by remote sensing (7, 12, 13, 23, 24). Excess gas is most pronounced in silicic magmas whose melts are characteristically S-poor, where the degree of excess degassing, defined as the ratio of measured SO2 yield to the petrologic method estimate, can range from 10 to 100 (11).

The “excess” S is predominantly supplied by an exsolved C-O-H-S vapor phase present within the magma before eruption (513). The important role of an S-rich vapor phase in determining total gas budgets of large silicic eruptions was demonstrated famously for the 1991 eruption of Pinatubo (25), with similar conclusions drawn for the Valley of Ten Thousand Smokes (8), Redoubt (9), Mount Saint Helens (26), and others. These findings are consistent with experimental studies on the fluid/melt partitioning of S in these characteristically S-poor rhyolitic melts (27).

Previous studies have established methods for quantifying a pre-eruptive magmatic gas phase at unmonitored eruptions (12, 16, 17) and have found that pre-eruptive gas can amount to up to 6 weight % (wt %) [~30 volume percent (volume %)] of crustal magma, implying that a significant proportion of erupted gas may be sourced from fluid generated during magma storage and evolution. Several lines of evidence, including the detection of substantial excess gas, suggest that silica-rich systems like Paektu may retain most of this vapor until eruption. Bubble migration through viscous rhyolitic melts is too slow for volatile transfer on eruptive time scales (13) and is inefficient particularly in crystal-poor melts (28). Magmatic convection, which is probably an important mechanism for pre-eruptive degassing of low-viscosity basaltic magmas (11), likely does not occur to an appreciable extent in kinetically sluggish high-viscosity silicic systems (13). Given these considerations, it is important to evaluate the presence, amount, and composition of any pre-eruptive fluid phase that may contribute to the total gas budget. The petrologic method, which only records degassing of dissolved volatiles and cannot account for the existence of a pre-eruptive vapor phase, may result in a significant underestimate of gas yields from large silicic eruptions, leading to inaccurate assessments of hazard and climate change potential.

Geologic background and the MI tephra

Paektu (42.0056°N, 128.0553°E) is an intraplate volcano whose 37 km2 summit crater is bisected by the political border between the Democratic People’s Republic of Korea (DPRK) and China. Despite its historical and geological significance, relatively little is known about Paektu, a volcano that has produced multiple large (volcanic explosivity index ≤ 6) explosive eruptions, including the ME, one of the largest volcanic events on Earth in the last 2000 years (Fig. 1).

The explosive ME deposited 23 ± 5 km3 dense rock equivalent (DRE) of material emplaced in two chemically distinct phases in the form of ash, pumice, and pyroclastic flow deposits (20): an initial nearly aphyric (≤3 volume % phenocrysts) comendite pumice (95% of tephra volume) and a late-stage phenocryst-rich (10 to 20 volume %) trachytic pumice (5% of tephra volume). The petrology and mineralogy of five ME pumices (four comendite and one trachyte; see Materials and Methods and table S1) plus major, trace, and volatile element concentrations in MGs and more than 100 comenditic and trachytic MIs hosted in anorthoclase, sanidine, clinopyroxene, olivine, and quartz were characterized by electron microprobe (EMP), Fourier transform infrared (FTIR) spectroscopy, and ion microprobe (SHRIMP; table S2). MIs used for this analysis are glassy and bubble-free, except for a handful of bubble-bearing MIs that were rehomogenized and analyzed for their CO2 contents (see Materials and Methods). None of the MIs used show signs of leakage or devitrification. Analyses of dissolved volatile concentrations reveal MI with moderate H2O (<4 wt %), minimal S (<300 ppm) and CO2 (<23 ppm in rehomogenized MI), and significant halogens (<5200 ppm Cl, <4000 ppm F).

To evaluate both comenditic and trachytic MI together, it is necessary to first establish that MI groups represent melts from the same magma lineage. To meet this criterion, the most evolved magmas (comendite) must represent residual melt generated by the fractional crystallization of the most primitive magmas (trachyte) as opposed to representing primitive magma contaminated by silicic country rock. The derivation of ME comendite from ME trachyte (sample CBS-TPUM; see Materials and Methods) was modeled by performing least squares linear regressions of the major element compositions of both glasses and natural crystal phases. The results of major element modeling with all major mineral phases are consistent with ME comendite being a ~21% residual melt (79% crystallization) from a CBS-TPUM parent (see Materials and Methods and tables S3 and S4). Modeled mineral abundances mimic those seen in natural pumice samples, with feldspar making up the large majority of the crystal population. The origin of Paektu comendites by fractional crystallization is also consistent with rare earth element concentrations and Pb, Nd, and Sr isotopic geochemistry (29), which suggests a lack of crustal contamination in causing the observed geochemical variations between rock types.

Trends in trace element concentrations in suites of MI can be used to confirm the extent of crystallization necessary to derive a daughter magma (at Paektu, comendite) from its parent (trachyte). During crystal growth, incompatible trace elements (for example, U, Cs, Nb, and Zr) will be excluded from fractionating crystal phases and become steadily enriched in the melt phase as magma evolution progresses. A perfectly incompatible element will increase in concentration in the melt as 1/(1 − X), where X is the weight fraction of crystals removed from the melt during fractional crystallization. Some elements will be less incompatible (for example, Zr in the presence of minor zircon). The element that shows the strongest enrichment is the best proxy for degree of crystal fractionation. The most incompatible element in the Paektu ME pumice suite is U, which increases from about 2.3 to 10.7 ppm from trachyte to most evolved comendite. Assuming crystal fractionation as the dominant driver of magma evolution, this 4.6-fold enrichment in U indicates ca. 78% crystallization from parental trachyte to evolved comendite melt, in agreement with major element regressions.


The incompatible element method for calculating volatile yields

With a parent-daughter relationship established for trachyte and comendite melts, each individual MI can be assumed to represent a discrete time step in the evolution of a single magma lineage. The presence of a pre-eruptive exsolved fluid phase can be indicated by examining trends in MI volatile contents throughout the lineage. Like incompatible elements, volatiles are largely excluded from growing crystal phases and become enriched in the melt as crystal fractionation progresses. Unlike other incompatibles, however, continued enrichment of volatiles can force fluid exsolution via “second boiling.” This occurs once the sum of the partial pressures of dissolved volatile species equals or surpasses local confining pressure. In the rock record, this process is recorded as the increase of volatiles and incompatible elements at a constant ratio until volatile solubility is reached, at which point the ratio between volatile and incompatible will begin to decrease. This incompatible element method for evaluating volatile exsolution in melts was pioneered by Anderson (30, 31) and refined by Wallace et al. (16, 17). Trends in volatiles versus U concentration in Paektu MI (Fig. 2 and Table 1) indicate fluid exsolution. The concentrations of H2O, S, and F in comendite MI (average of 2.40 wt %, 110 ppm, and 3354 ppm) are much lower relative to their expected concentrations after crystal fractionation in a fluid-absent system (6.18 wt %, 905 ppm, and 4220 ppm), indicating the partitioning of these elements out of the melt.

Fig. 2 Incompatible element (U) versus volatile plots illustrating trends in Paektu MI chemistry.

An increase in U concentration represents continued crystallization-evolution of the magma. The “Fluid absent crystallization” line represents the expected trend for MI chemistry in a fluid-absent system, where volatile and incompatible elements will be equally enriched in the melt as crystallization progresses. The red lines (in essence, the difference between the expected and actual average comendite MI volatile concentration) represent the amount of that volatile exsolved as a separate fluid (bubble) phase. The orange dashed lines (in essence, the difference between average comendite MI and MG concentrations) represent the amount of that volatile degassed during decompression and ascent upon eruption. Error bars are calculated on the basis of relative errors reported in table S2.

Table 1 Volatile concentrations in MI and MG, expected fluid-absent concentrations in evolved MI, and Δvolatile values for pre- and syn-eruptive fluids.

U concentration is used as an index of differentiation. All analyses are in parts per million, unless noted. Numbers in parentheses represent 1σ uncertainty in units of the last reported digits. H2O and CO2 are measured by transmission or ATR FTIR. S, F, and U are measured by SHRIMP. Cl is measured by EMP. n, number of MI averaged.

View this table:

The difference between expected and actual volatile concentrations in the most evolved melts is a function of the amount of the volatile species lost from the melt. We estimate that 3.78 wt % H2O, 795 ppm S, and 896 ppm F were lost from the melt during differentiation. We can then calculate the absolute masses of H2O, S, and F lost from the melt asEmbedded Image(1)where M is the volatile mass in Tg, Δvolatile is the weight % of the volatile lost from melt, ρ is the density of the melt (2.4 × 1012 kg/km3), Φ is the crystal-free fraction of the magma (0.97), and V is the erupted magma volume in km3 DRE (22.8-km3 comendite). Assuming that all volatiles lost from the melt partitioned into a fluid phase, this gives maximum pre-eruptive fluid masses of H2O, S, and F of approximately 2006, 42, and 46 Tg, respectively.

Because ME pumice contains both S- and F-bearing crystal phases (sulfide, molybdenite, fluorapatite, and fluorite), the volatile masses we calculate here are actually partitioned between fluid and solid, and the actual fluid mass will depend critically on the extent to which solid phases sequester S and F. In addition, S and F may be added to the fluid budget if volatile-bearing crystal phases formed in the trachytes (in essence, before MI formation) were broken down during eruption, a process inferred at other systems (32). Calculations of solid sequestration and breakdown are discussed in detail in the “Contributions from solid phases” section.

Fluid phase modeling of Cl and CO2

The relative masses of the remaining species Cl and CO2 in the fluid phase may be calculated on the basis of their concentrations in MI using experimentally established partitioning behaviors if key magma chamber parameters (P and T) are known.

The observation that volatile concentrations in comendite MI fall below expected concentrations along a fluid-absent crystallization path (Fig. 2) suggests that ME comendites were fluid-saturated at the time of MI entrapment. H2O solubility models can be used to indicate minimum saturation pressures of comendite MI, which, in a fluid-saturated system, are indicative of MI trapping pressures. MI H2O contents range from ~1 to 4 wt %, indicative of pressures from ~150 to 900 bars [using the solubility model of Liu et al. (33)], with most of the MI clustering between 200 and 600 bars (average of ca. 400 bars; table S4). This corresponds to a magma chamber depth between 0.5 and 3.5 km, assuming an average crustal density of 2.4 g/cm3. Note that the maximum magma chamber depth of 3.5 km is consistent with seismicity recorded during the 2002–2005 period of seismic unrest at Paektu (34).

Comendite storage temperature can be estimated using the results of experimental phase equilibria determinations for volatile-bearing Paektu rocks. The natural phase assemblage of ME comendites was reproduced at ≤1000 bars and below ~720°C. This relatively low temperature is consistent with titanium-in-zircon and alkali feldspar-glass geothermometry in ME comendites (~740 ± 40°C) (35) and experimental studies on similar peralkaline silicic rocks (36). The low crystal contents of comendite pumices indicate near-liquidus conditions or a temperature near 720°C at 400 bars.

At low pressures (≤1300 bars), chloride brines exsolved during second boiling will immediately separate into a Cl-rich hydrosaline liquid and a Cl-poor vapor phase (37). This makes the calculation of the amount of Cl in any vapor exsolved from a melt less straightforward, because the method used for H2O, S, and F (Eq. 1) does not account for any Cl taken up in a hydrosaline liquid and would thus overestimate the amount of Cl in the vapor. The composition of coexisting liquid brine and vapor phases in the H2O-NaCl system can be determined using the phase relations established at magmatically relevant P and T (Fig. 3) (38). Degassing of Cl during the ME was relatively minimal; the vapor phase confined to <2 wt % NaCl (1.2 wt % Cl), with any remaining exsolved Cl taken up by the hydrosaline liquid.

Fig. 3 Isothermal P-X projection showing the compositions of coexisting vapor and liquid brine phases in the H2O-NaCl system after the study by Bodnar et al. (38).

The solvus that exists at multiple temperatures below ca. 2000 bars indicates the composition of both Cl-poor vapor (left of the critical curve) and Cl-rich liquid brine (right of the critical curve), which coexist as immiscible fluids at magmatic pressures and temperatures. The yellow star indicates the NaCl content of pre-eruptive ME fluid at P and T inferred for comenditic ME melt inclusions, and the blue bar indicates the range in the NaCl content of the vapor given the range in P calculated on the basis of MI H2O contents.

The concentration of CO2 in the fluid phase can be calculated using mixed volatile solubility models appropriate for ME rocks, given the known concentrations of H2O and CO2 dissolved in the melt. Pristine (bubble-free) ME MIs contain no detectable CO2 down to a detection limit of ±2 ppm via FTIR. Because CO2 will often exsolve within MIs to form a CO2-rich vapor bubble (39), a population of bubble-bearing MIs was chosen for rehomogenization, in which MI-bearing crystals were held at magmatic P and T for a time sufficient to dissolve the vapor bubble back into the melt (also see Materials and Methods). Rehomogenized MIs were then measured by FTIR and were found to contain up to 23 ppm (re)dissolved CO2.

Because of its low solubility, even relatively low concentrations of CO2 in the melt can have a large effect on the equilibrium fluid-phase composition. VolatileCalc (40) and MagmaSat (41) predict that a fluid in equilibrium with an ME comendite melt containing 23 ppm CO2 and ~2.5 wt % H2O would contain ~11 wt % CO2 at P around 400 bars, assuming an ideal mixing behavior between all species.

With the absolute masses of H2O, S, and F in the fluid known and the relative masses of Cl and CO2 known, we can now solve for the total fluid mass asEmbedded Image(2)

This gives a total pre-eruptive fluid mass of 2386 Tg, made up of 84.1 wt % H2O (2006 Tg), 1.8 wt % S (42 Tg), 1.9 wt % F (46 Tg), 0.8 wt % Cl (20 Tg), and 11.4 wt % CO2 (272 Tg; Fig. 4 and Table 2).

Fig. 4 Breakdown of fluid source and composition.

This diagram illustrates the proportional contribution and composition of both pre- and syn-eruptive fluid, plus the total gas yield, which is equal to the sum of the two fluid types. Colored circles represent each fluid species (H2O, S, F, Cl, and CO2), and the area of the circles corresponds to the fluid mass in teragrams (see figure legend). Gray bars illustrate the proportional contribution of each fluid type (pre- and syn-eruptive) to the total yield and correspond to the vertical axis.

Table 2 Pre- and syn-eruptive and total gas compositions and masses.
View this table:

Syn-eruptive degassing

During ascent and decompression of an erupting magma, additional gas will be released from the melt caused by the drop in solubility of volatiles (particularly H2O) as confining pressure decreases. The total gas yield of the ME is thus the sum of the pre-eruptive fluid and syn-eruptive degassing. The latter can be estimated using the so-called petrologic method as the difference between pre-eruptive volatile contents as recorded in MIs and in degassed MGs, which, unlike MIs, are not enclosed within rigid crystals and thus respond to changes in pressure during ascent. Syn-eruptive fluid masses can be calculated by modifying Eq. 1, where F is the syn-eruptive fluid mass in Tg and Δvolatile is the difference between average MI and MG volatile concentrations in weight %. Values for the ME are given in Table 2 and shown in Fig. 4, which illustrates the composition of both fluid types, the proportional contribution of each fluid type to the total eruptive yield, and the proportion of each volatile species contributed by each fluid type.

Although most of the erupted gas was H2O (84.1 wt %, 3121 Tg), with up to 45 Tg of S and 78 Tg of halogens (58 Tg of F and 20 Tg of Cl), the ME is among the largest contributors of climate-forcing volcanic gas to the atmosphere in the last 2 ka (Fig. 1). Note that more than 90% of S released during the ME was supplied by a pre-eruptive fluid phase (42 Tg) and very little from syn-eruptive degassing (2.7 Tg), consistent with findings of previous studies that suggest that S emissions from silicic magmas are largely sourced from pre-eruptive fluid (811, 2527). Our syn-eruptive S estimate of 2.7 is similar to an earlier S yield calculated for the ME using the petrologic method (2 Tg) (20), a technique known for underestimating total volatile yields because of its inability to consider contributions from pre-eruptive fluid.

Contributions from solid phases

S- and F-bearing minerals molybdenite (42), Cu-Fe-sulfide, fluorite, and apatite have been identified in Paektu rocks, and so, any volatiles taken up by these phases need to be taken into account when calculating the composition of the fluid phase. Molybdenite (MoS2) is rarely found in contact with glass and is most commonly found included in feldspar crystals. Most feldspars are inclusion-free, but those that do contain molybdenite inclusions often contain many (1020). Cu-Fe-sulfide is extremely rare and typically found included in mafic phenocrysts. Fluorite is extremely rare and only appears as small anhedral crystals with disequilibrium textures included in feldspar and pyroxene. Fluorapatite is less rare and is found in contact with glass and included in feldspar and pyroxene. The lack of Ba-rich rims on feldspars or Ti-rich rims on quartz is inconsistent with magmatic defrosting, suggesting that crystals in ME pumices were grown and stored at high temperatures and were not derived from remelted crystal mush (43, 44).

S-bearing minerals are rare in ME pumice (probably, most of the crystals were separated into an unerupted crystal mush), and so, the total mass of precipitated sulfide and molybdenite is unconstrained. However, we can estimate the total mass of each of these phases precipitated during fractionation from trachyte to comendite by examining trends in their constituent elements (Cu in the case of sulfide and Mo for molybdenite). Given ~79% fractionation of dominantly feldspar to produce ME comendite from trachyte, Cu and Mo in the melt, if perfectly incompatible, will be expected to increase in concentration by a factor of 4.6 during fractionation. For example, a concentration of 8 ppm Cu in bulk trachyte results in an expected concentration of 36.8 ppm Cu in comendite MI. Actual concentrations in comendite MI average 14.5 ppm. Using the approach as in Eq. 1, we calculate that 1.22 Tg of Cu was removed from the melt during fractionation. Assuming that all of this went into Cu-Fe-sulfide (composed of 35 wt % Cu, 32.5 wt % Fe, 32.5 wt % S), this equates to a total sulfide mass of 3.49 Tg (~0.007 wt %) and a total sequestration by sulfide of 1.1 Tg of S.

A similar approach can be taken for molybdenite. Given concentrations of 5 to 9 ppm Mo in bulk trachyte result in an expected concentration of 30.6 ppm Mo in comendite. We did not measure Mo in MIs, but we can use the Mo concentration measured in comendite bulk rock via x-ray fluorescence (XRF) as a proxy because most of the Mo will be contained in the melt. Bulk comendite Mo concentrations of 10 ppm give us a Mo shortfall of 20.6 ppm. Assuming that all of this went into molybdenite with a concentration of 60 wt % Mo and 40 wt % S, this equates to a total molybdenite mass of 1.82 Tg (~0.004 wt %) and a sequestration by molybdenite of 0.7 Tg of S. This would lower our pre-eruptive S to 40.2 Tg and our total S yield to 43.2 Tg. Because Cu and Mo are both fluid mobile, some of the Cu and Mo lost from the melt probably went into a vapor phase. Thus, these calculations represent the maximum possible S sequestered by S-bearing crystals.

The above approach is not readily applicable to F-bearing crystals because their constituent elements partition heavily into silicate phases. For example, in the case of fluorite, the necessary assumption that 100% of any Ca shortfall in comendites was sequestered by fluorite is invalid in the presence of Ca-bearing silicates. Instead, we assume the precipitation of 0.01 wt % fluorite and 0.05 wt % fluorapatite. Given an average concentration of 48.7 wt % F in fluorite and 4 wt % F in fluorapatite, this equates to 0.00487 and 0.002 wt % F sequestered in the two phases, respectively. Again, using the same approach as in Eq. 1, this gives 3.6 Tg of F sequestered in solid phases, lowering our total pre-eruptive F to 42.4 Tg and our total F yield to 54.4 Tg.

The possibility for sulfides to volatilize during eruption, releasing their S into the fluid, as has been inferred at other systems, is worth noting (32). Any S release from the syn-eruptive breakdown of solid phases that precipitated during evolution from trachyte to comendite will be accounted for by the calculation of S release via Eq. 1, and so, we assume this contribution to be negligible.

Uncertainties in gas yield calculations

Uncertainties have been estimated for each element in chemical analyses and all additional parameters used to calculate volatile fluxes (see Materials and Methods). A Gaussian propagation of uncertainties through all equations (see the Supplementary Materials) results in relative 1σ uncertainties in H2O, S, F, Cl, and CO2 total mass yields of ±37%, 22%, 22%, 45%, and 41%, respectively, with a 41% error on the total gas yield. Errors are affected most strongly by the large uncertainty in estimated erupted magma volume (±21%). The large errors on Cl and CO2 are due to the calculation of their pre-eruptive mass by mass balance and include propagated error associated with the pre-eruptive masses of H2O (±27%), S (±22%), and F (±22%).

In addition to uncertainty on the accuracy of our calculations, our reported fluxes can be evaluated by comparing fluid and melt volatile concentrations to experimentally established fluid-melt partition coefficients at relevant P, T, and fo2. The compositions of coexisting quartz, fayalitic olivine, and ilmenite in natural comendite pumices indicate storage at an oxygen fugacity just below that of the quartz-fayalite-magnetite (QFM) buffer. Using the program QUILF (45), we obtain a logfo2 of QFM-1.2 at 720°C and 400 bars. The partitioning behavior of S between fluid and melt (DSfl/m) has been experimentally determined for Paektu-like comendites (ME comendite MI average Na2O + K2O/Al2O3 = 1.22) from the Greater Olkaria Volcanic Complex at 800°C and 1500 bars: samples ND and SMN with (Na2O + K2O)/Al2O3 = 1.05 and 1.31, respectively (46). At an fo2 near QFM-1, Scaillet and Macdonald report a sulfur DS of 250 ± 125 for sample SMN, the most Paektu-like composition. Our estimates of 1.8 wt % S in the fluid and 110 ppm S in the melt yield a DS of 163, within the 50% error on the experimental value. An additional 8 Tg of S fluid from solid-phase breakdown gives a total pre-eruptive S yield of 50 Tg, bringing the DS to 191.


All possible gas: A maximum gas yield for the ME

Our method for calculating volcanic gas yields considers all possible gas generated during magma storage and evolution. This may represent a protracted time span (zircon ages and Ra/Th isotopes suggest a comendite residence time of ca. 6 to 9 ky) (47, 48), and so, it is likely that some of the pre-eruptive gas left the system via passive degassing before the ME. This much lower energy degassing would not inject volatiles into the stratosphere and so would not play any significant role in climate forcing. Our reported volatile masses thus represent maximum possible gas yields from the ME. Previous estimates, calculated using the petrologic method (20), represent minimum possible yields, and so, the true gas yield of the ME must lie somewhere between these values.

Over the past few decades, there has been mounting evidence to suggest that substantial pre-eruptive fluid is common in large silicic systems and that much of this fluid is retained until eruption (812). The very detection of excess S at these systems, whereby the petrologic method underestimates true gas yields by one to two orders of magnitude (11), demonstrates the ability of silicic magmas to retain very large volumes of gas during storage. This has been explained by the inefficiency of volatile diffusion and bubble transport through kinetically sluggish, high-viscosity melts (13). More recently, results from laboratory experiments have demonstrated that, although vapor bubbles may migrate freely through crystal-rich magmas, they tend to accumulate in crystal-poor magmas because of the interplay between capillary stresses and the viscosity contrast between melt and vapor (28).

If substantial pre-eruptive fluid was retained until the ME, then our maximum yields provide a more accurate foundation for evaluating potential climate impacts. Although our data provide compelling evidence for the formation of a pre-eruptive vapor phase, without direct evidence of the retention of this vapor, it may be just as likely that the petrologic estimate is a more appropriate tool in the case of the ME. Because pre-eruptive vapor is interpreted to play such a substantial role in many explosive silicic eruptions, however, the methods we present here may serve as a way forward in evaluating the full range of possible gas yields from silicic eruptions that have previously been evaluated using only the petrologic method.

Internal triggering of explosive eruptions by gas overpressure

It is commonly interpreted that evolved crystal-poor rhyolites represent segregated regions of melt surrounded by crystal-rich mush (~50 volume % crystals) and outer rigid sponge (>65 volume % crystals) zones (Fig. 5) (49). Accumulation of the residual melt is likely driven by crystal settling and compaction and/or filter pressing of a crystal mush [for example, studies of Bachmann and Bergantz (50), Bea (51), and Dufek and Bachmann (52)]. This is consistent with the extensive fractionation (ca. 79 wt % for Paektu comendites) required to generate such silica-rich melts.

Fig. 5 Schematic drawing of the ME magma chamber beneath Paektu before eruption.

An upper chamber of crystal-poor comendite magma is generated by crystal settling and/or filter pressing of a crystallizing comendite magma (50, 51, 75). Crystals settle out of the main chamber and condense into mush (~50% crystals) and outer rigid sponge (>65% crystals) zones (47). Parental trachyte feeder dikes supply the comendite chamber with new magma, promoting prolonged crystallization-driven gas exsolution within the chamber. Exsolution of volatile-bearing melt during ascent further propels the magma upward and, along with pre-eruptive gas, results in a ca. 25-km-high Plinian eruption column and extensive tephra fallout. S (45 Tg) and halogens (78 Tg) are released, and much of this gas is likely injected into the upper atmosphere.

Assuming that all the pre-eruptive fluid is retained until eruption, 2386 Tg of fluid accounts for about 6 wt % of the comendite magma. This value is similar to pre-eruptive fluid fractions of 1 to 6 wt % determined for other intermediate to silicic eruptions of >0.01 km3 of magma, such as Pinatubo, the upper Bishop Tuff, Rabaul, El Chichon, and others (12, 16, 17). The high viscosity of ME comendite magma may have stifled pre-eruptive gas loss by hampering the formation of interconnected bubble networks before eruption and probably promoted vertical gradation in the mass fraction of gas such that the most gas-rich magma resides in the roof of the reservoir (17). At Paektu, the total calculated gas yield equates to ~30 volume % of the magma, near the threshold above which bubble networks become sufficiently interconnected to allow permeable gas flow (53). The buildup of a gas-rich cap within the magma chamber due to prolonged fluid exsolution could create the conditions for explosive fragmentation of the magma (54), promoting triggering of the ME by internal overpressure. Gas exsolution triggering of the ME is consistent with explosive eruption dynamics inferred from tephra, which initiated with a 25-km-high Plinian eruption column (20) recorded as widespread and voluminous fall deposits of highly vesiculated (~75 volume %) pumice clasts.

The top-down emptying of the comendite chamber is evidenced by the small-volume late-stage trachyte, likely representative of a mafic feeder dike (possibly one of many) responsible for delivering parental melt into a growing comendite reservoir. At Paektu, the lack of chemically distinct rims on crystals indicates insufficient time for an injection of trachyte magma to mix into the comendite reservoir and disturb chemical equilibria. This suggests that the eruption was not triggered by the energy from a single mafic injection but rather by the buildup of overpressure due to prolonged input, differentiation, and exsolution of parental magma.

Climatic impact of the ME

Taken alone, the significant S and halogen yields from the ME suggest that it had the potential to affect Northern Hemisphere atmosphere and surface temperatures in the years following the eruption. Climatic impacts of ancient eruptions can be difficult to quantify and depend on many factors, including the volcano’s latitude and season of eruption, in addition to total S yield. For example, perturbations in climate following eruptions at high latitudes are much smaller and shorter-lived than if the same eruption had occurred in the tropics (55). This effect is compounded if the eruption occurred during winter, when the removal of stratospheric aerosols in polar regions is enhanced (56). Thus, this suggests that climate forcing after the ME, which occurred at 42°N latitude and likely during Northern Hemisphere winter (based on tree-ring estimates) (57), may have been diminished.

Signatures of climate forcing after large eruptions can sometimes be identified as sharp increases in sulfate deposited in polar ice cores, which preserve a timeline of paleoclimate chemistry going back several thousand years (58, 59). Historical records detailing possible climatic effects after the ME (and records of the eruption itself) are ambiguous (60, 61), and so, sulfate loading as determined from ice core records is currently the best proxy for determining climatic impact. A spike in sulfate deposited in Greenland ice around 940 CE has been positively linked to the ME by matching chemical fingerprints in simultaneously emplaced tephra to ME ash (19). The relatively small sulfate load (9 kg km−2) is concluded by Sun et al. (19) to be indicative of minimal climate forcing from the ME when compared to sulfate deposited after Tambora (~40 kg km−2) and Krakatau (~15 to 18 kg km−2), two tropical eruptions that were followed by drops in global surface temperatures (18, 62). This conclusion is shared by other authors (56, 60), all of whom cite low S emissions, based on the 2-Tg estimate of Horn and Schmincke (20), as an important factor along with seasonality and high latitude of the Paektu eruption.

Our new gas yield estimate considers all possible gas emitted during the ME (all pre-eruptive vapor is assumed to be retained until eruption), and although this represents a maximum value, we stress that the range of total possible volatile yields must be considered in silicic magma systems, where the exsolution and retention of a pre-eruptive vapor phase may be substantial. Given our new upper bound for the S yield from the ME, we suggest that if large-scale surface cooling after the eruption was minimal, it was likely due to seasonal and latitude controls on the transport of sulfate aerosols to the arctic, as proposed by Kravitz and Robock (56), and is not reflective of S release from Paektu.


Major element fractional crystallization modeling

Our analysis of the volatile contents in trachytic and comenditic MIs requires that both MI groups be derived from the same parent magma and that little to no contamination occurred during the derivation of a comendite liquid from a trachyte parent. The relationship between rock types and the weight fraction of residual comendite liquid from a trachytic parent was estimated by linear least-squares multiple regression, mass-balancing the major oxide composition of the assumed trachyte parent with the compositions of comendite glass plus its phenocryst phases (63). The compositions used for this modeling are given in table S3.

The regression was designed to model the relationship between dependent (parent) and independent (daughter liquid plus crystals) variables by fitting a linear equation to observed data. The general solution is expressed as Data = Fit + Residual, where the Fit term gives phase proportions for independent variables and the Residual term represents the variance between the calculated parent liquid composition and the actual (observed) parent liquid. A model run was considered successful when (i) the sum of all phase proportions was close to 1, (ii) all phase proportions were positive, and (iii) the residuals were small (≪1).

For this modeling, we used the Microsoft Excel plug-in StatPlus and a python script using the pandas and statsmodels.api libraries (available at Both of these methods gave the same results. The two most successful model runs are given in table S4. Both predicted trachyte glass compositions, in good agreement with the observed composition, and both were consistent with the derivation of PEK-62 comendite as a 21% residual melt from a CBS-TPUM parent.

Samples and sample preparation

We analyzed MGs, phenocrysts, and more than 100 glassy MIs hosted in feldspar, pyroxene, quartz, and olivine from four ME comendite pumices collected from airfall deposits in the DPRK (PKTU-TP, PEK-26, PEK-56, and PEK-62) and one ME trachyte pumice collected from an airfall deposit in China (CBS-TPUM; table S1). All samples are registered in the SESAR database with International GeoSample Numbers (IGSNs) IACPKTUTP, IAC00001L, IAC00002D, and IAC00002J for comendites and IGSN IAC00000H for the trachyte.

Comendite pumices represent the initial volumetrically dominant (95 volume % of tephra) phase of the eruption. Comendites are nearly aphyric with ≤3 volume % phenocrysts of mainly anorthoclase feldspar, hedenbergite, fayalitic olivine, ilmenite, and rare quartz plus accessory fluorapatite, chevkinite, molybdenite, Cu-Fe-sulfide (intermediate solid solution), and zircon. All comendite pumices used in this study were sampled from airfall deposits.

Representative pumice clasts were taken from all levels of the ME comendite pumice deposits around the Paektu volcano. Deposits were typically 1 to 4 m thick and clast-supported without interstitial ash. Clasts with diameters of 1 to 10 cm were crushed, and phenocrysts were picked by hand. Care was taken to only choose crystals containing MIs that appeared to be pristine. Any MI not fully enclosed by the host phenocryst, inclusions that showed evidence of post-entrapment crystallization, inclusions touching cracks in the host phenocrysts, or inclusions with bubbles were rejected. A population of bubble-bearing (but otherwise pristine) inclusions was chosen for rehomogenization experiments and subsequent CO2 concentration analysis.

Trachyte CBS-TPUM represents the final stage of eruption and makes up only 5 volume % (vol %) of the entire tephra volume (20). CBS-TPUM was sampled from an airfall pumice deposit in China (42.137474°N and 128.350082°E) and belongs to the late-stage trachytic eruptive member, which is composed of phenocryst-rich (10 to 20 vol %) trachytic pumice emplaced as pumice fall (outcrops mainly in China) and proximal agglutinates (mantling crater walls) followed by slightly welded trachytic pyroclastic flows (mainly in DPRK). The phenocryst assemblage in the trachyte consists of sanidine, hedenbergite, ilmenite, and rare olivine plus accessory fluorapatite and chevkinite. Aggregates of xenocrystic feldspar, quartz, zircon, and magnetite are also present. MI-bearing phenocrysts were selected from CBS-TPUM in the same manner as for comendite pumices. Volatiles in MIs were measured using FTIR (H2O and CO2), EMP (Cl), and SHRIMP (F and S).

Fourier transform infrared spectroscopy

Transmission FTIR measurements were performed at the U.S. Geological Survey (USGS) (Menlo Park) using a Thermo Scientific Nicolet iN10 MX mapping FTIR. MI-bearing crystals were doubly polished so that the MI was exposed on both sides, resulting in crystals with thicknesses ranging from ~30 to 150 μm. Transmission infrared spectra were obtained in the wave number range of 6000 to 1000 cm−1 using a KBr beam splitter with a variable aperture size between ~25 and 150 μm2. The sample stage was continuously purged with N2, and a background measurement was taken before each spectrum, allowing for determinations in CO2 content down to a detection limit of ca. 2 ppm.

H2O was retrieved by measuring the peak height of either the total water peak at 3550 cm−1 or the combination of the OH and molecular water peaks at 4500 and 5200 cm−1, respectively. CO2 was retrieved by measuring the peak height at 2350 cm−1. The following absorption coefficients for appropriate rhyolitic melts were chosen from the literature: 1.73 [OH (64)], 1.61 [molecular H2O (64)], 80 [total water (65)], and 945 [CO2 (65)]. Calculations of H2O and CO2 concentrations from FTIR data and associated errors are given in the Supplementary Materials.

Comendite MG and several pyroxene crystals were singly polished and analyzed for H2O using attenuated total reflectance (ATR) FTIR (66). The technique was calibrated by measuring the peak heights at 3450 cm−1 for several rhyolite standard glasses from Lowenstern and Pitcher (66) with varying known water concentrations. Water concentrations in unknowns were then calculated using the calibration curve constructed from standards. ATR measurements are accurate to within about 0.2 wt %.

Electron microprobe

Major element and Cl concentrations in MI, MG, and crystal phases were analyzed by EMP at the USGS (Menlo Park) using a JEOL 8900 Superprobe. Glasses were measured using an accelerating voltage of 15 kV and a defocused beam of either 5 or 10 μm. Na was measured at a beam current of 2 nA, a peak counting time of 10 s, and a background counting time of 5 s. Major and minor elements were measured with a beam current of either 8 or 20 nA, respectively, a peak counting time of 20 s, and a background counting time of 10 s.

Significant Na loss during EMP measurement caused by the migration of Na+ ions away from the electron beam is a common problem associated with microprobe analyses of highly alkaline, hydrous glasses [(67) and references therein]. Na loss was confirmed in hydrous ME comendites by comparing Na contents in crystal-free hydrous glasses synthesized from ME comendite PEK-62 pumice with those by XRF whole-rock analysis. Na loss was found to scale linearly with H2O content of the melt asEmbedded Image(3)where “Apparent Na2O loss” is the difference between Na2O in the whole rock and that measured with EMP (normalized to 100% on an anhydrous basis). A correction for Na2O using this equation was applied to all MIs with known H2O contents.

Sensitive high-resolution ion microprobe

Trace elements, F, and S in MI and MGs were measured at the USGS-Stanford SHRIMP-RG facility at Stanford University following the analysis and data reduction procedures as outlined by Wright et al. (68). Standards used were Macusani glass, NIST glasses 611 and 613 (69), and RLS glasses 132, 140, 158, and 37 (70).

Volatile concentrations and CO2 determination

Comendite MIs have moderate H2O contents (≤4 wt %), low S (110 ppm), and substantial halogen concentrations (Cl ≤ 4000 ppm, F ≤ 3500 ppm). No measureable CO2 was present in pristine comendite or trachyte MI. During post-entrapment cooling of MI-bearing crystals, the relatively large thermal contraction of melt compared to crystal host can force exsolution of a vapor bubble (shrinkage bubble) within the MI (71). These shrinkage bubbles can contain large amounts of CO2 because of its relatively low solubility. A subset of bubble-bearing comendite MI was used in high-P rehomogenization experiments to check whether bubbles contained substantial CO2. No trachytic MIs were rehomogenized. Several crystals were placed in a titanium-zirconium-molybdenum–type pressure vessel at 1000 bars and 900°C for approximately 20 min and then rapidly quenched. Postrun MIs were bubble-free (in essence, the bubbles had dissolved back into the melt) and showed a small increase in CO2 (≤23 ppm). The lack of abundant dissolved CO2 is consistent with low MI trapping pressures (average of ca. 400 bars) inferred based on comendite H2O contents. At these low pressures, a silicic magma will be expected to have degassed most, if not all, of its CO2, which exsolves from silicate melt much deeper than other magmatic volatile species.


Average uncertainties were determined for values obtained via EMP, FTIR, and SHRIMP and are reported as relative fractional error in table S2. Reported errors for major oxides and Cl by EMP represent the average SD for all MIs. For each oxide, this was determined by calculating the SD on multiple analyses for a given MI (each MI has n number of analyses; table S2) and then averaging the SDs for all MIs and dividing by the average value for a given oxide.

Uncertainties in trace elements, F, and S by SHRIMP were determined on the basis of the average SD of values obtained during multiple measurements of standards. Uncertainties in H2O and CO2 by FTIR were determined by propagating errors associated with the measurement of sample thickness, spectral peak heights, and calculation of sample density. All values, associated errors, and final propagated error for each MI are reported in the Supplementary Materials.


Supplementary material for this article is available at

table S1. Whole-rock pumice compositions by XRF, microprobe, and FTIR.

table S2. Major, trace, and volatile element analyses in ME MIs.

table S3. Compositions of phases used in least squares linear regression modeling.

table S4. Results of least squares linear regression modeling of the derivation of PEK-62 comendite from CBS-TPUM trachyte.

table S5. H2O contents in comenditic MI and modeled saturation pressures.

table S6. Excel file listing raw FTIR data, all values necessary for calculation of H2O and CO2 concentrations from FTIR data, and associated errors of measured and calculated values.

table S7. Excel file listing all values used to calculate volatile fluxes and associated errors.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: This work is a part of the international Mt. Paektu Geoscientific Group project. It received strong support from the Pyongyang International Information Centre of New Technology and Economy (PIINTEC), the American Association for the Advancement of Science (AAAS), the Royal Society, and the Environmental Education Media Project (EEMP). We especially thank K. M. Ryol (PIINTEC); N. Neureiter and R. Stone (AAAS); M. Poliakoff, L. Clarke, C. Dynes, and B. Koppelman (Royal Society); and J. Liu (EEMP) for their logistical support and backing. K.I. thanks J. Vasquez and M. Coble for assistance in performing SHRIMP measurements and analysis, L. Hayden for microprobe and scanning electron microscopy analysis, and S. Burgess for help in shaping the manuscript. All authors thank two anonymous reviewers and the editorial support of B. Schoene. Funding: K.I. was supported by the NSF under award no. 1349486 and by AAAS. Fieldwork was supported by the Richard Lounsbery Foundation. Author contributions: K.I., C.O., and J.O.S.H. conceived the research project. K.I. and A.D. performed analytical work. T.S. and J.L. supported K.I.’s position at the USGS and contributed substantially to the content of the manuscript and to data collection. K.I., K.J.-S., R.K.-H., J.J.-N., S.K.-H., H.S.-H., C.O., J.O.S.H., A.D., K.W.L., and R.K.-R. were integral in sample collection. All authors assisted in the preparation of the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.
View Abstract

Stay Connected to Science Advances

Navigate This Article