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

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. INTRODUCTION 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 weatheringinduced 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; KrissansenTotten 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 Paleo zoic 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 1GSA Data Repository item 2019219, detailed methods, supplementary figures and data tables, and associated references, is available online at http:// www .geosociety .org /datarepository /2019/, or on request from editing@ geosociety .org. CITATION: Soreghan, G.S., Soreghan, M.J., and Heavens, N.G., 2019, Explosive volcanism as a key driver of the late Paleozoic ice age: Geology, v. 47, p. 600–604, https:// doi .org /10 .1130 /G46349.1 Manuscript received 19 January 2019 Revised manuscript received 5 April 2019 Manuscript accepted 9 April 2019 https://doi.org/10.1130/G46349.1 © 2019 The Authors. Gold Open Access: This paper is published under the terms of the CC-BY license. Published online 2 May 2019 Downloaded from https://pubs.geoscienceworld.org/gsa/geology/article-pdf/47/7/600/4775200/600.pdf by guest on 06 July 2019 Geological Society of America | GEOLOGY | Volume 47 | Number 7 | www.gsapubs.org 601 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. GLACIATION AND VOLCANISM IN THE LPIA 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 vol canic 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 stagelevel 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 200 220 240 260 280 300 320 340 360 380 400 N um be r

Although pCO 2 was relatively low during the LPIA, it is difficult to ascribe the LPIA to pCO 2 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 2 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 2 levels inferred for the late Paleozoic (Royer, 2014).Also, mismatches persist between reconstructed pCO 2 and the glacial record.For example, pCO 2 (Foster et al., 2017) exhibits an ambiguous relationship to the timing of onset, demise, and peak of the LPIA: whereas pCO 2 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 2 lows ca.305 Ma and 298 Ma, closer to peak ice conditions, but depict pCO 2 rising at the apex of the LPIA (ca.298-295 Ma; Fig. 1).Moreover, climate and climate-ice sheet models indicate a CO 2 glaciation threshold at ~560 ppmv (Lowry et al., 2014), but high-resolution pCO 2 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 2 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 2 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 Paleo zoic glacioeustasy" (Horton and Poulsen, 2009, p. 715), in which order-of-magnitude variations in pCO 2 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 2 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 2 and indirect forcings related to volcanism.

GLACIATION AND VOLCANISM IN THE LPIA
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 vol canic 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

RADIATIVE CONSEQUENCES OF VOLCANIC AEROSOLS
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 5 yr.In contrast, the main Quaternary SLIP (Taupo Volcanic Zone, New Zealand; Bryan, 2007) erupts through tively 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 vol canics in similar areas and intervals are identified as ignimbrites associated with cal deras (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 2 .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 2 sulfate that accounts for the additional sulfate radiative forcing at the reconstructed pCO 2 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 2 emission since 1979 CE (Schmidt et al., 2018).
Our results indicate that volcanic sulfate forcing maintains effective pCO 2 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 2 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 2 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 2 (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 longerperiod maxima in pCO 2 , 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 2 and formation of larger, less-reflective sulfate aerosol particles (e.g., Timmreck et al., 2010).

INDIRECT EARTH-SYSTEM EFFECTS OF SULFATE AEROSOLS
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 2 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 2 , 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.

Figure 2 .
Figure 2. Spatial and temporal distribution of volcanism through study interval (see methods in the Data Repository for details, and Table DR2 for source data [see footnote 1]).A: Scatter plot of ages of volcanic deposits (400-200 Ma) versus reconstructed paleolatitudes, designated by type of volcanic deposit.B: Two-dimensional histogram of age of vol canic deposits and reconstructed paleolatitudes for each point, presented as a heat map.Histogram values are smoothed by a 5 × 5 gaussian filter.Peak intensity occurs at ca. 300-290 Ma, with greatest contributions occurring in equatorial and (secondarily) southern mid-latitudes.
Figure 3. Radiative forcing from volcanic sulfate and CO 2 for 400-200 Ma.A: Reconstructed pCO 2 (black line; Foster et al., 2017) versus equivalent pCO 2 (red dotted line) accounting for sulfate radiative forcing (RF); horizontal blue line shows model pCO 2 threshold (560 ppm) for glaciation (Lowry et al., 2014).Gray bar is interval for which records of glaciation exist (Fig. 1; Table DR1 [see footnote 1]).LPIAlate Paleozoic ice age.B: Comparison of negative radiative forcing attributable to pCO 2 (blue line), negative radiative forcing attributable to volcanic sulfate aerosol (red line), and total radiative forcing (yellow line) relative to 560 ppmv CO 2 .Histogram at base shows frequency of glacial deposits (Fig. 1).