Atmospheric CO2 exerts a robust and well-documented control on Earth’s climate, but the timing of glaciation during the late Paleozoic Ice Age (LPIA; ca. 360–260 Ma) is inconsistent with pCO2 reconstructions, hinting at another factor. Stratospheric volcanic aerosols produce a large but temporary negative radiative forcing under modern conditions. Here we examine explosive volcanism over 200 m.y. of Earth history to show that the LPIA corresponded with a sustained increase in volcanism in both tropical and extratropical latitudes. A major peak in explosive volcanism at ca. 300 Ma likely corresponded to stratospheric sulfur-injecting eruptions at least three to eight times more frequent than at present. This level of volcanism created a steady, negative radiative forcing potentially sufficient to initiate and, most critically, sustain icehouse conditions, even under increasing levels of pCO2, and helps resolve discrepancies between glacial timing and CO2 records. Accounting for the radiative forcing effects of CO2 and sulfate indicates that both are required to explain the LPIA, with sulfate producing an especially strong effect at peak icehouse ca. 298–295 Ma. Frequent explosive volcanism would have increased atmospheric acidity, enhancing the reactivity of iron in abundant volcanic ash and glacially generated mineral dust, thus strengthening the climate impact of volcanism through a marine biological pump further primed by feedback with glaciation.

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 pCO2 (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 CO2 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 Repository1). 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 pCO2 along with lower solar luminosity (Crowley and Baum, 1992) is the preferred explanation for the LPIA (Montañez et al., 2016), although causes for CO2 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 CO2 (e.g., Feulner, 2017) and accelerated silicate weathering. Plant expansion and CO2 drawdown, however, are asynchronous, leading others to propose weathering-induced drawdown related to Pangean orogeny (Goddéris et al., 2017).

Although pCO2 was relatively low during the LPIA, it is difficult to ascribe the LPIA to pCO2 alone. Both enhanced silicate weathering and organic carbon burial are negative feedbacks (Walker et al., 1981; Berner, 2004; Krissansen-Totten and Catling, 2017). Lower pCO2 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 pO2 levels inferred for the late Paleozoic (Royer, 2014). Also, mismatches persist between reconstructed pCO2 and the glacial record. For example, pCO2 (Foster et al., 2017) exhibits an ambiguous relationship to the timing of onset, demise, and peak of the LPIA: whereas pCO2 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 pCO2 lows ca. 305 Ma and 298 Ma, closer to peak ice conditions, but depict pCO2 rising at the apex of the LPIA (ca. 298–295 Ma; Fig. 1). Moreover, climate and climate–ice sheet models indicate a CO2 glaciation threshold at ∼560 ppmv (Lowry et al., 2014), but high-resolution pCO2 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 pCO2 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 pCO2 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 pCO2 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 pCO2 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 pCO2 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 = log10[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 >105 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 pCO2. 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 pCO2sulfate that accounts for the additional sulfate radiative forcing at the reconstructed pCO2 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 SO2 emission since 1979 CE (Schmidt et al., 2018).

Our results indicate that volcanic sulfate forcing maintains effective pCO2 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); pCO2 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 CO2 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 pCO2 (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 pCO2, 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 SO2 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 pCO2 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 pCO2, 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.

1GSA Data Repository item 2019219, detailed methods, supplementary figures and data tables, and associated references, is available online at, or on request from
Gold Open Access: This paper is published under the terms of the CC-BY license.