Explosive volcanism as a key driver of the late Paleozoic ice age

Earth’s climate has fluctuated between icehouse and greenhouse states, characterized by the presence or absence of continental ice and strongly correlated with variations in pCO₂ (Berner, 2004; Royer et al., 2004; Jagoutz et al., 2016; McKenzie et al., 2016). The late Paleozoic Ice Age (LPIA; ca. 360–260 Ma) archives the longest Phanerozoic icehouse, and its sedimentary record preserves evidence for atmospheric CO₂ intermittently as low as that of Quaternary glaciations (Montañez et al., 2007, 2016). Glaciation began in the Late Devonian (Isaacson et al., 2008; Lakin et al., 2016); after an Early–Middle Mississippian minimum, glaciation peaked in the earliest Permian (Fig. 1A; Table DR1 in the GSA Data Repository¹). At times, continental ice reached latitudes as low as 32° (Evans, 2003), and alpine glaciation is hypothesized for equatorial uplands (Soreghan et al., 2014). The LPIA ultimately collapsed—Earth’s only example of icehouse termination on a fully vegetated planet (Gastaldo et al., 1996).

Diminishing pCO₂ along with lower solar luminosity (Crowley and Baum, 1992) is the preferred explanation for the LPIA (Montañez et al., 2016), although causes for CO₂ drawdown are disputed. Long-standing research (Berner, 2004; Royer et al., 2004) posits the primacy of the terrestrial carbon cycle, wherein evolution of plants sequestered CO₂ (e.g., Feulner, 2017) and accelerated silicate weathering. Plant expansion and CO₂ drawdown, however, are asynchronous, leading others to propose weathering-induced drawdown related to Pangean orogeny (Goddéris et al., 2017).

Although pCO₂ was relatively low during the LPIA, it is difficult to ascribe the LPIA to pCO₂ alone. Both enhanced silicate weathering and organic carbon burial are negative feedbacks (Walker et al., 1981; Berner, 2004; Krissansen-Totten and Catling, 2017). Lower pCO₂ can weaken the first effect by reducing precipitation and temperature, and the second by reducing net primary productivity through increased plant respiration (Pagani et al., 2009; Gerhart and Ward, 2010), particularly at the higher pO₂ levels inferred for the late Paleozoic (Royer, 2014). Also, mismatches persist between reconstructed pCO₂ and the glacial record. For example, pCO₂ (Foster et al., 2017) exhibits an ambiguous relationship to the timing of onset, demise, and peak of the LPIA: whereas pCO₂ reached a nadir at ca. 338–334 Ma, peak glaciation occurred ca. 298–295 Ma (Fig. 1). The highest-resolution reconstructions (Montañez et al., 2016) cover a brief interval of the LPIA (ca. 311–298 Ma) and show pCO₂ lows ca. 305 Ma and 298 Ma, closer to peak ice conditions, but depict pCO₂ rising at the apex of the LPIA (ca. 298–295 Ma; Fig. 1). Moreover, climate and climate–ice sheet models indicate a CO₂ glaciation threshold at ∼560 ppmv (Lowry et al., 2014), but high-resolution pCO₂ reconstructions for the interval near peak icehouse (Montañez et al., 2016) show a high-frequency oscillation both above and below this threshold. Finally, climate models cannot account for hypothesized equatorial glaciation (Soreghan et al., 2014) in moderate-elevation uplands without invoking pCO₂ levels (<200 ppmv) that would stress modern vegetation (Pagani et al., 2009), calling into question either the data or the modeling.

Decoupling between glaciation and pCO₂ could reflect a hysteresis effect in which higher positive radiative forcing from insolation and/or greenhouse gases is required to ablate an ice sheet than the negative radiative forcing necessary to form it (Zhuang et al., 2014). However, such a hysteresis creates a “paradox of late Paleozoic glacioeustasy” (Horton and Poulsen, 2009, p. 715), in which order-of-magnitude variations in pCO₂ in addition to orbital forcing are required to generate the glacioeustasy of the LPIA. But an analogous paradox for the Miocene Antarctic ice sheet is now ascribed to oversimplification of ice-sheet geophysics (Gasson et al., 2016). Indeed, the ice-sheet hysteresis effect was reduced for late Paleozoic simulations by introduction of dynamic vegetation (Horton and Poulsen, 2009). In short, such hysteresis appears to be an emergent property of some models, and thus fails to explain the shortcomings in correlation of the glaciation and pCO₂ records.

Here we integrate geologic data and radiative calculations to explore the hypothesis that the onset, acme, and, especially, prolonged extent of the LPIA were driven in part by unusually intense explosive volcanism prevalent during Pangean assembly, operating in concert with pCO₂ and indirect forcings related to volcanism.

To assess the role of volcanic aerosols, we compiled data (see the Data Repository) on (1) glacial deposits (Table DR1) and (2) radioisotopically dated felsic to intermediate volcanic units globally from 400 to 200 Ma (Fig. 1; Table DR2), enveloping the LPIA and bounding greenhouse intervals. We record one glacial deposit per sedimentary basin (cf. Hambrey and Harland, 1981) to minimize bias, and use stage-level binning to reflect the dating resolution of most deposits. Glacial deposits are reported for every stage of the time scale from the Famennian to the Wuchiapingian (Fig. 1). For the volcanic compilation, we focus on felsic to intermediate compositions, as these are the main source of eruptions that inject sulfur (S) species into the stratosphere (Mather and Pyle, 2015). Once in the stratosphere, S species form sulfate aerosols that linger over multiple years and reflect solar radiation, producing a negative radiative forcing (Robock, 2000; Crowley, 2000). In addition, volcanic eruptions can impact climate on longer time scales by cooling the ocean and affecting the initial thermal state of intermediate and deep water masses; that is, the ocean integrates the sub-decadal radiative perturbations of multiple, individual eruptions into a sustained forcing on centennial to millennial time scales (e.g., Stenchikov et al., 2009; Zanchettin et al., 2012; Miller et al., 2012). Extreme stratospheric S injection also may have occurred in some continental flood basalt eruptions (Glaze et al., 2017), but these eruptions have recurrence time scales similar to deep ocean overturning scales (∼2000 yr) (Kasbohm and Schoene, 2018) and so would produce a less-sustained, although high-magnitude, forcing (Schmidt et al., 2016).

Frequency of volcanic units varies by more than an order of magnitude (Fig. 1). This variability remains evident even when considering only ignimbritic units, typically associated with proximal deposition of magnitude (M = log₁₀[erupted volume, in kg] − 7) ≥6.5 eruptions (Engwell and Eychenne, 2016). Greenhouse times of the Triassic and pre–Late Devonian exhibit very low frequencies of explosive volcanism, as does the Early–Middle Mississippian (ca. 359–330 Ma)—a time with significantly lower frequency of glacial deposits. The first substantial increase in volcanism corresponds with the Late Devonian glaciation ca. 360 Ma. Acceleration in volcanism corresponds with the Late Mississippian onset of large-scale Gondwanan glaciation (ca. 330 Ma), and peak glaciation corresponds to peak volcanism, which is strongly concentrated along the equatorial Central Pangean Mountains and the southern mid-latitude Kennedy-Connors-Auburn silicic large igneous province (SLIP) of Australia (Fig. 2).

Peak volcanism likely corresponds to a level of activity (frequency of explosive eruptions at all magnitudes) at least three to eight times greater than that of the last 2.5 m.y. (Fig. 1D). This is calculated by ratioing the number of ignimbrites to the expected value if the Quaternary record of ≥M6.5 eruptions were sampled by our compilation (see the Data Repository). Scaling volcanic activity to ignimbrites is justifiable if large silicic eruptions primarily occurred within SLIPs. Extrapolating frequency of larger eruptions from the frequency of smaller Holocene eruptions might underestimate the number of eruptions ≥M8 (Deligne et al., 2010), which suggests that M8 eruption frequency is statistically and perhaps dynamically decoupled from most eruptions that inject significant S into the stratosphere. This decoupling might arise because the largest (almost M9) eruptions require formation of deep magma chambers over >10⁵ yr. In contrast, the main Quaternary SLIP (Taupo Volcanic Zone, New Zealand; Bryan, 2007) erupts through relatively young crust, shifting its eruptive history toward smaller volumes and shorter recurrence, as expected of a more continuous size-frequency eruptive distribution (Wilson et al., 1984).

Uncertainties in volcanic activity relative to that at present (Fig. 1D) include (1) extent of ignimbrite deposits from the volcanic vent and (2) whether we could resolve distinct eruptions over the same 2.5 m.y. interval. For example, our estimate assumes (1) that large, shallow eruptions spread ignimbrites farther than large, deep eruptions or small eruptions, and (2) resolution of only the largest eruption of an individual volcano within a 2.5 m.y. interval. Because of these last two assumptions, our estimate is a lower (conservative) bound. M8 eruptions are thought to be mostly identified and preserved in the Quaternary, but the majority of M7 eruptions probably are not (Brown et al., 2014), suggesting a gradual erasure of the record, particularly of smaller eruptions.

This range of estimated volcanic activity relative to the present day still may be too conservative. Many Pennsylvanian–early Permian volcanic rocks in Western Europe (paleo-equator) dated according to our strict standards are generically classified as pyroclastics (see detailed methods in the Data Repository) (Wilson et al., 2004; Paulick and Breitkreuz, 2005), but volcanics in similar areas and intervals are identified as ignimbrites associated with calderas (Geißler et al., 2008; Willcock et al., 2013). Uncertainties in the completeness and characterization of the stratigraphic record are addressed in the detailed methods in the Data Repository.

Recent work using detrital zircons and arc length as proxies for Neoproterozoic–Phanerozoic continental arc volcanism (McKenzie et al., 2016; Cao et al., 2017) inferred a lull in arc volcanism during the LPIA and linked this to reduced pCO₂. Our findings do not refute this, and indeed are reinforced by the raw detrital zircon data of McKenzie et al. (2016), which depict the Carboniferous as having the maximal peak in zircon yield of the Phanerozoic. Because our data derive from volcanics dated at high precision, they incorporate no lag time (normally associated with exhumation of arc granitoids), and thus integrate both arc and other (e.g., transform, intraplate) volcanism. Also, our data strongly suggest that explosive (ignimbritic) volcanism predominated during peak LPIA (ca. 300 Ma), associated with orogenic collapse in the Variscan-Hercynian belt (Wilson et al., 2004) of the paleo-equator and the convergent-to-transform transition in the Kennedy-Connors-Auburn SLIP of the southern mid-latitudes (Murray et al., 1987) (Fig. 2).

A volcanic radiative forcing was estimated and also converted to a pCO₂^sulfate that accounts for the additional sulfate radiative forcing at the reconstructed pCO₂ level (Fig. 3). The key assumptions underlying this calculation are that (1) average volcanic activity from 500 BCE to 1900 CE is similar to that over the last 2.5 m.y. (Brown et al., 2014; Toohey and Sigl, 2017), (2) volcanic sulfate aerosol optical depth (SAOD) averages 0.010 over that period (Toohey and Sigl, 2017), and (3) the transfer function between SAOD and radiative forcing is linear and equivalent to an efficiency of –23.7 W m^–2, which agrees well with most estimates based on latest-generation climate modeling and satellite-based estimates of volcanic SO₂ emission since 1979 CE (Schmidt et al., 2018).

Our results indicate that volcanic sulfate forcing maintains effective pCO₂ near or below the 560 ppmv glaciation threshold from ca. 360 to 257 Ma, corresponding exactly with the interval recording glacial deposits of the LPIA (Fig. 3A); pCO₂ did not go below 560 ppmv until ca. 345 Ma, and rose above it as early as ca. 290 Ma, well before icehouse termination. Although CO₂ forcing likely drove the cold of the Late Mississippian, we posit that volcanic forcing was particularly critical for sustaining cold conditions thereafter (Fig. 3). The interval of maximum frequency of glacial deposits (ca. 298–295 Ma) falls in the midst of maximum (negative) radiative forcing (–1.6 W m^–2) induced by sulfate aerosols (Fig. 3). Indeed, this level of negative radiative forcing might have compensated for the positive radiative forcing relative to 560 ppmv (0.4–3.2 W m^–2) due to high-frequency fluctuations in pCO₂ (600–920 ppmv) during the Late Pennsylvanian (Gzhelian; 304–299 Ma) reported by Montañez et al. (2016). Thus, by acting as a contrary radiative forcing to Milankovitch-scale as well as longer-period maxima in pCO₂, sulfate aerosols could have sustained icehouse conditions from late Carboniferous well into Permian time.

Note that the transfer function between SAOD and radiative forcing responds to two factors that may have been different during the LPIA relative to modern: (1) planetary albedo and (2) the oxidizing power of the atmosphere. Higher planetary albedo (e.g., bright, bare soil or greater cloudiness) would reduce the net reflective effect of sulfate aerosol. Greater oxidizing power would enable faster oxidation of SO₂ and formation of larger, less-reflective sulfate aerosol particles (e.g., Timmreck et al., 2010).

Direct radiative forcing from explosive volcanism was potentially high enough to have stabilized the LPIA, but also highly uncertain in magnitude. However, a significant increase in explosive volcanism also could have enabled higher carbon burial rates, helping to maintain low pCO₂ levels. First, volcanic ash is a critical source of nutrients like Si, Fe, and P to otherwise nutrient-limited marine and terrestrial ecosystems (Lee et al., 2018). Second, the attendant atmospheric acidity would have enhanced iron solubility of both ash and mineral dust (Oakes et al., 2012). This may explain the high concentrations of highly reactive iron in atmospheric dust of the late Paleozoic, which would have enhanced marine primary productivity and thus carbon burial (Sur et al., 2015) and perhaps fostered the exceptional prevalence of red beds and associated sulfates in this interval.

As we move forward on Earth with rapidly increasing concentrations of atmospheric pCO₂, stratospheric aerosol geoengineering is increasingly cited as a means to mitigate climate change (Wigley, 2006). But lessons from deep time can shed light on the spectrum of outcomes, that what begins with a reduction of sunlight reaching the surface may not end there. In this sense, the LPIA—as Earth’s only example of icehouse collapse on a vegetated planet—is particularly intriguing as a case study.

G. Soreghan acknowledges funding from National Science Foundation (NSF) Sedimentary Geology and Paleobiology (SGP) grant EAR-1338331. M. Soreghan acknowledges funding from NSF SGP grant EAR-1053018. N. Heavens acknowledges funding from NSF SGP grant EAR-1337463. We thank X. Qi, K. Baczkowski, A. Sweet, and A. Oordt for help with data archival; I. Montanez for sharing in-press age data; and reviewers (C.-T. Lee, W. DiMichele, and three anonymous) and editor J. Parrish for constructive comments on previous versions of this manuscript.