Influence of temporally varying weatherability on CO₂-climate coupling and ecosystem change in the late Paleozoic

Earth's penultimate icehouse period, the late Paleozoic ice age (LPIA), was a time of dynamic glaciation and repeated ecosystem perturbation, which was under conditions of substantial variability in atmospheric pCO₂ and O₂. Improved constraints on the evolution of atmospheric pCO₂ and O₂∕CO₂ ratios during the LPIA and its subsequent demise to permanent greenhouse conditions are crucial for better understanding the nature of linkages between atmospheric composition, climate, and ecosystem perturbation during this time. We present a new and age-recalibrated pCO₂ reconstruction for a 40 Myr interval (∼313 to 273 Ma) of the late Paleozoic that (1) confirms a previously hypothesized strong CO₂–glaciation linkage, (2) documents synchroneity between major pCO₂ and O₂∕CO₂ changes and compositional turnovers in terrestrial and marine ecosystems, (3) lends support for a modeled progressive decrease in the CO₂ threshold for initiation of continental ice sheets during the LPIA, and (4) indicates a likely role of CO₂ and O₂∕CO₂ thresholds in floral ecologic turnovers. Modeling of the relative role of CO₂ sinks and sources active during the LPIA and its demise on steady-state pCO₂ using an intermediate-complexity climate–carbon cycle model (GEOCLIM) and comparison to the new multi-proxy CO₂ record provides new insight into the relative influences of the uplift of the Central Pangean Mountains, intensifying aridification, and increasing mafic rock to granite rock ratio of outcropping rocks on the global efficiency of CO₂ consumption and secular change in steady-state pCO₂ through the late Paleozoic.

1 Introduction

Earth's penultimate and longest-lived icehouse period (340 to 290 Ma) occurred under the lowest atmospheric CO₂ concentrations of the last half-billion years (Foster et al., 2017) and, potentially, the highest atmospheric pO₂ of the Phanerozoic (Glasspool et al., 2015; Krause et al., 2018; Lenton et al., 2018). Anomalous atmospheric composition, along with 3 % lower solar luminosity (Crowley and Baum, 1992), may have primed the planet for a near-miss global glaciation (Feulner, 2017). Notably, Earth's earliest tropical forests assembled and expanded during this icehouse period (the late Paleozoic ice age; LPIA), leading to the emergence of large-scale wildfire. Paleotropical terrestrial ecosystems underwent repeated turnovers in composition and architecture, culminating in the collapse of wetland (coal) forests throughout tropical Pangea at the close of the Carboniferous (Cleal and Thomas, 2005; DiMichele, 2014), possibly promoting the diversification and ultimate dominance of amniotes (Pardo et al., 2019). In the marine realm, global rates of macroevolution (origination, extinction) decreased, in particular among tropical marine invertebrates, and genera with narrow latitudinal ranges went extinct at the onset of the LPIA (Stanley, 2016; Balseiro and Powell, 2019; Powell, 2005). Low marine macroevolutionary rates continued through to the demise of the LPIA in the early Permian (Stanley and Powell, 2003; McGhee, 2018).

Reconstructions of late Paleozoic atmospheric pCO₂ document a broad synchroneity between shifts in CO₂, glaciation history, glacioeustasy, and restructuring of paleotropical biomes, underpinning the hypothesized greenhouse gas forcing of sub-million-year glacial–interglacial cycles (Montañez et al., 2016) and the terminal demise of the LPIA (Montañez et al., 2007). For late Paleozoic pCO₂ (and pO₂) reconstructions, however, broad intervals of low temporal resolution and significant uncertainties limit the degree to which mechanistic linkages between atmospheric composition, climate, and ecosystem change can be further evaluated. Moreover, the potential impact of large-magnitude fluctuations in atmospheric O₂∕CO₂, which characterized the late Paleozoic, on the biosphere has been minimally addressed. On longer timescales (≥10⁶ years), the relative roles of potential CO₂ sinks and sources on secular changes in late Paleozoic atmospheric CO₂, and in turn as drivers of glaciation and its demise, remain debated (McKenzie et al., 2016; Goddéris et al., 2017; Macdonald et al., 2019).

Here, we present a multiproxy atmospheric pCO₂ reconstruction for a 40 Myr interval (313 to 273 Ma) of the late Paleozoic, developed using new leaf-fossil-based estimates integrated with recently published age-recalibrated Pennsylvanian pCO₂ estimates of 10⁵-year resolution (Montañez et al., 2016) and re-evaluated fossil-soil-based (paleosol) CO₂ estimates for the early Permian (Montañez et al., 2007). Our new multiproxy record offers higher temporal resolution than existing archives while minimizing and integrating both temporal and CO₂ uncertainties. This pCO₂ reconstruction, together with new O₂∕CO₂ estimates of similar temporal resolution, permits refined interrogation of the potential links between fluctuations in atmospheric composition, climate shifts, and ecosystem events through Earth's penultimate icehouse period. Moreover, comparison of the new 40 Myr CO₂ record with modeled steady-state pCO₂ and seawater ⁸⁷Sr∕⁸⁶Sr over the same interval provides new insight into the relative importance and evolution of CO₂ sinks and sources during late Paleozoic glaciation and its turnover to a permanent greenhouse state.

2 Materials and methods

A brief account of the methods is presented here; more details are presented in the Supplement. Primary data generated or used in this study are deposited in the Dryad Digital Repository (Richey et al., 2020) and can be accessed at https://doi.org/10.25338/B8S90Q.

2.1 Sample collection and analysis

To build the pCO₂ record, 15 plant cuticle fossil species/morphotypes were used, which were collected from eight localities in Illinois, Indiana, Kansas, and Texas, USA, including four well-studied Pennsylvanian interglacial floras (Sub-Minshall (313 Ma; Šimůnek, 2018), Kinney Brick (305.7 Ma; DiMichele et al., 2013), Lake Sarah Limestone (303.7 Ma; Šimůnek, 2018), and Hamilton Quarry (302.7 Ma; Hernandez-Castillo et al., 2009a, b, c); Figs. 1a and S2–4 in the Supplement; Richey et al., 2020; https://doi.org/10.25338/B8S90Q). The Pennsylvanian estimates were integrated into a previously published pCO₂ reconstruction (313 to 296 Ma; Montañez et al., 2016) of 10⁵-year resolution built using pedogenic carbonates and wet-adapted seed fern fossils (Figs. 2b, S1b). The Permian estimates were integrated with previously published latest Carboniferous and early Permian pedogenic carbonate-based CO₂ estimates (Montañez et al., 2007), derived from paleosols from successions throughout Arizona, New Mexico, Oklahoma, Texas, and Utah, USA (Fig. 1a; Richey et al., 2020; https://doi.org/10.25338/B8S90Q). The pedogenic carbonates and leaf cuticle fossils span a broad region of Pennsylvanian and early Permian tropical Euramerica (Fig. 1b). Ages of samples used in Montañez et al. (2007, 2016) were recalibrated and assigned uncertainties using the latest geologic timescale (Ogg et al., 2016) and biostratigraphic and geochronologic controls (see the Supplement; Richey et al., 2020; https://doi.org/10.25338/B8S90Q).

Figure 1 Sampling localities in present-day and late Paleozoic geographic context. (a) Sampling locations of pedogenic carbonates and plant fossils and their position relative to the late Pennsylvanian (310 and 300 Ma) and early Permian (290 to 270 Ma) equator (the colors of the flora localities correspond to that of the paleo-equator at that time). White band traversing NM (New Mexico) and CO (Colorado) is the area of inferred shortening during the Laramide and Sevier orogenies. Map modified from Montañez et al. (2007). (b) Earliest Permian (290 Ma) paleogeography (Scotese, 2016); shading corresponds to paleotopographic/bathymetric scale on the right. Inset box is the location of panel (a).

Figure 2 Late Paleozoic CO₂ estimates. (a) New and revised (Montañez et al., 2007) pCO₂ estimates, bootstrapped LOESS trend, and 75 % confidence interval (CI). Revised pedogenic carbonate-based estimates were made using δ¹³C_OOM (blue filled circles; n=28; Fig. S1) and δ¹³C_POM (open black circles; n=16; Fig. S1). Trend line is the median of 1000 bootstrapped LOESS analyses; dashed intervals indicate low data density and higher uncertainty. See Fig. S1 for error bars on individual CO₂ estimates and the 95 % CI, and see Richey et al. (2020) (http://doi.org/10.25338/B8S90Q) for the full dataset. (b) Multiproxy CO₂ record and individual estimates (this study and age-recalibrated values by Montañez et al., 2016; n=165), documented glacial deposits (Soreghan et al., 2019), occurrence of peats (coal), and best estimate of timing (and uncertainties) of magmatic episodes: 1a = Tarim 1, China (∼300 Ma); 1b = Tarim 2 (292–287, peak ∼290 Ma); 1c = Tarim 3 (284–272, peak ∼280 Ma; Chen and Xu, 2019); 2 = Skagerrak-centered, NW Europe (297.5±3.8 Ma; Torsvik et al., 2008); 3a = Panjal Traps, NW India (289±3 Ma; Shellnutt, 2018); 3b = Qiangtang Traps, Tibet (283±2 Ma; Zhai et al., 2013); 4 = Choiyoi, W Argentina (beginning 286.5±2.3 Ma, continuing for up to 39 Myr; Sato et al., 2015). Trend lines are as in (a); dashed intervals across the Carboniferous–Permian boundary (298.9 Ma) indicate overlap of the two LOESS trend lines.

Plant cuticle and organic matter occluded within pedogenic carbonates (OOM) were rinsed in 3 M HCl to remove carbonates and analyzed at the Stable Isotope Facility, University of California, Davis, using a PDZ Europa ANCA-GSL elemental analyzer interfaced to a PDZ Europa 20-20 IRMS (isotope-ratio mass spectrometer). External precision, based on repeated analysis of standards and replicates, is $<\pm \mathrm{0.2}$ ‰. For Hamilton Quarry (HQ), all material was previously mounted on slides for taxonomic analysis (Hernandez-Castillo et al., 2009a, b, c). Because of this, biomarker δ¹³C values of bulk stratigraphic sediment samples were used (Richey et al., unpublished data; see the Supplement). HQ n-C₂₇₋₃₁ n-alkane δ¹³C was analyzed using a Thermo Scientific GC Isolink gas chromatograph connected to a Thermo Scientific MAT 253 mass spectrometer. Standard deviation of n-alkane δ¹³C was ±0.3 ‰. For biomarker δ¹³C, a +4 ‰ correction was used to account for fractionation during biosynthesis (Diefendorf et al., 2015), and the standard deviation of all values was used as the uncertainty (1.6 ‰, which is five times the analytical precision).

2.2 Models

The MATLAB model Paleosol Barometer Uncertainty Quantification (PBUQ; Breecker, 2013), which fully propagates uncertainty in all input parameters, was used to derive pedogenic carbonate-based CO₂ estimates (Figs. 2a, S1a). For each locality, paleosols of inferred different soil orders were modeled separately. We applied improved soil-specific values for soil-respired CO₂ concentrations (S_(z); Montañez, 2013) and the δ¹³C of organic matter occluded within carbonate nodules (δ¹³C_OOM; Fig. S5) as a proxy of soil-respired CO₂ δ¹³C. For samples where OOM was not recovered, estimates were revised using PBUQ and the plant fossil organic matter δ¹³C used in Montañez et al. (2007) (δ¹³C_POM; Fig. S5). Because of the limited amount of carbonate nodules remaining after study by Montañez et al. (2007), δ¹³C_OOM was substituted for δ¹³C_POM for localities that occur in the same geologic formation, and a large error (±2 ‰) was used to account for the uncertainty in this approach. PBUQ model runs conducted in this study resulted in a small subpopulation of biologically untenable CO₂ estimates (i.e., ≤170 ppm; Gerhart and Ward, 2010). To limit estimates below that threshold, two changes to the PBUQ MATLAB code were made (see the Supplement for details). All other input parameters remained unchanged from Montañez et al. (2007).

For cuticle fossil-based (Figs. S2–4) CO₂ estimates (Figs. 2a, S1a), we utilized a mechanistic (non-taxon-specific) gas-exchange model (Franks et al., 2014). For some fossil cuticles, pore length (PL) was measured directly; for others, PL was inferred from guard cell length (GCL; Table S1 in the Supplement). Guard cell width was estimated via GCL using the prescribed gymnosperms and ferns scaler (0.6; Franks et al., 2014; Table S1).

For both stomatal- and pedogenic carbonate-based CO₂ modeling, we calculated δ¹³C of atmospheric CO₂ using the carbonate δ¹³C record generated from an open-water carbonate slope succession (Naqing succession, southern China; Buggisch et al., 2011), contemporaneous estimates of mean annual temperature (Tabor and Montañez 2005; Tabor et al., 2013), and temperature-sensitive fractionation between low-Mg calcite and atmospheric CO₂ (Romanek et al., 1992; Eq. S2; Table S1).

We used the spatially resolved intermediate-complexity GEOCLIM model (Goddéris et al., 2014) to quantitatively evaluate how steady-state atmospheric CO₂ may have responded to changes in weatherability and the relative influences of different CO₂ sources and sinks. The spatial distributions of the mean annual runoff and surface temperature were calculated offline for five time increments (Goddéris et al., 2017), covering the period of interest, and for various atmospheric CO₂ levels using the 3-D ocean–atmosphere climate model FOAM and the approach as described in Donnadieu et al. (2006). GEOCLIM uses generated lookup tables to calculate steady-state atmospheric CO₂ for a given continental configuration and to account for paleogeography and relief. Although the GEOCLIM model does not include an explicit surface distribution of lithology, the weathering rate of mafic rocks and continental granites are calculated using different methods, and the impact of physical erosion on granite weathering is accounted for (Goddéris et al., 2017). For mafic surfaces, a simple parametric law is used, linking the surface of the considered grid cell, local runoff, and mean annual temperature to the local mafic weathering rate. The calibration of the GEOCLIM model was performed at the continental-scale by tuning the parameters of the model so that 30 % of the alkalinity generated by the weathering of silicates originates from the weathering of mafic rocks (GEOCLIM REG; Dessert et al., 2001; Goddéris et al., 2014).

2.3 O₂∕CO₂ ratio

O₂∕CO₂ ratios (Fig. 3a) were calculated using the 10 000 CO₂ estimates produced by our modeling and combined with O₂ estimates obtained using geochemical mass balance and biogeochemical models (Krause et al., 2018; Lenton et al., 2018). Unreasonably high O₂∕CO₂ ratios (generally those that correspond to CO₂≲200 ppm) were removed from the resulting O₂∕CO₂ dataset.

Figure 3 Late Paleozoic O₂∕CO₂ ratios and pCO₂, as well as a comparison to environmental and biotic events. (a) O₂∕CO₂ estimates made using CO₂ values of this study and averaged time-equivalent modeled O₂ (Krause et al., 2018; Lenton et al., 2018). Trend line is the median of 1000 bootstrapped LOESS analyses; gray horizontal line is present-day O₂∕CO₂. (b) Bootstrapped Pennsylvanian and Permian LOESS analyses (from Fig. 2a), with significant overlap across the Pennsylvanian–Permian boundary interval. The shading indicates CO₂ above (orange) and below (shades of blue) the mean value for the 16×10⁶-year record through the late Pennsylvanian reported in Montañez et al. (2016). Temporal changes in terrestrial (c) and marine (d) ecosystems. Plant biomes from Montañez (2016): wetland biome (WB) 2 (i.e., cordaitalean/lycopsid co-dominance), WB 3 (i.e., lycopsid dominated with cordaitaleans decreasing and tree ferns increasing in occurance and prominence), WB 4 (i.e., tree fern dominated), dryland biome (DB) 1 (i.e., cordaitalean dominated), DB 2 (i.e., walchian dominated). Diagonal arrows indicate 10⁵-year glacial–interglacial shifts between wet- and dry-adapted floras.

2.4 Statistical analyses

We utilize a bootstrap approach that assesses uncertainties of both CO₂ (or O₂∕CO₂) and age. Each age uncertainty was truncated to ensure no overlap in locality ages, which were constrained by their relative stratigraphic position to one another (Richey et al., 2020; https://doi.org/10.25338/B8S90Q). The 10 000 modeled CO₂ estimates were trimmed by 28 % to remove anomalously high or low values. The means of the resulting 7200 CO₂ estimates were compared to the trimmed means of the 10 000 CO₂ estimates to ensure that trimming did not alter the central tendency of the data. Locality ages were resampled and perturbed, assuming that the individual ages and truncated age uncertainties represent the mean and standard deviation of the ages. Similarly, the trimmed CO₂–O₂∕CO₂ datasets were resampled with replacement, and the resampled ages and estimates were used to build 1000 resampled datasets. Each resampled dataset was subjected to LOESS analysis (0.25 smoothing), and the median and 95 % and 75 % confidence intervals of the LOESS analyses were calculated (Figs. 2, 3a, b, S1). The Pennsylvanian and Permian portions of the record were analyzed separately due to differing data density, with significant overlap across the Pennsylvanian–Permian boundary interval (Figs. 2b, 3b, S1b).

Figure 4 Analysis of the statistical significance of short-term CO₂ fluctuations. (a) White intervals (A–K) delineate short-term highs or lows in the CO₂ LOESS trend used for binning (n=11; bins ±0.5 to 1 Myr resolution). Raw stomatal-based and pedogenic carbonate-based CO₂ estimates generated by Monte Carlo analysis (10 000 model runs per CO₂ estimate; data in shaded intervals were not used). CO₂ between bins was compared by calculating the mean of the lowest through to the 10 000th (highest) Monte Carlo values for all CO₂ points in each bin and comparing the means of the two bins sequentially. (b–h) Histograms of the percent change between each of the 10 000 Monte Carlo means of the adjacent bins. Negative values indicate a decrease in value between bins and positive values an increase. The number above each histogram bar is the “percent change” values represented in each bar. The percent of the 10 000 model runs that confirm a given increase or decrease in the LOESS trend is indicated by the % value shown on the right side of each panel. See Sect. 2 for further details.

To test the validity of short-term fluctuations in the LOESS CO₂ trend, we undertook further analysis of the raw Monte Carlo data produced by PBUQ and the mechanistic stomatal model in several short-term increments. Eleven short-term highs or lows (A–K in Fig. 4a) were designated and used to form bins of ±0.5 to ±1 Myr. Within an individual bin, each shown “bootstrapped” CO₂ data point is the trimmed mean of 10 000 Monte Carlo model runs. The Monte Carlo model runs for each data point were sorted from lowest to highest CO₂ value, and the lowest CO₂ values for each data point within the bin were averaged. This averaging was repeated sequentially for each of the 10 000 values creating 10 000 means for each bin (n=11). To evaluate whether a visually perceived rise or fall (e.g., A to B decrease or B to C increase) is statistically valid, the 10 000 means of two adjacent bins were compared sequentially with one another (i.e., the mean of the lowest values of one bin was compared to the mean of the lowest values of the adjacent bin) in order to calculate a percent change $\left(\right((V_{\mathrm{2}}-V_{\mathrm{1}})/V_{\mathrm{1}})\times \mathrm{100})$ for each of the 10 000 model runs, resulting in 10 000 percent changes for each set of adjacent bins. The percent of the 10 000 comparisons that confirm an increase or decrease between bins is reported (Fig. 4b–j) as a measure of the statistical significance of the short-term fluctuations in CO₂ concentration visually observed on the LOESS trend.

3 Results

Revised early Permian mineral-based CO₂ estimates define a substantially narrower range (45–1150 ppm; Fig. 2a) than previous estimates (175–3500 ppm) made using the same pedogenic carbonate sample set (Montañez et al., 2007) while maintaining the original trends and including fewer photosynthetically untenable concentrations (≤170 ppm; Gerhart and Ward, 2010). New early Permian cuticle-based estimates show a high level of congruence by locality and broad plant functional type, falling within the revised pedogenic carbonate-based CO₂ range (Figs. 2a, S1a). Similarly, stomatal-based estimates for the four Pennsylvanian interglacial floras are within the estimated pCO₂ range defined by the pedogenic carbonates (Figs. 2a, S1a) and late-glacial wetland plant fossils (Montañez et al., 2016). Notably, the newly integrated record confirms elevated atmospheric CO₂ concentrations (482 to 713 ppm ($-\mathrm{28}/+\mathrm{72}$ ppm)) during Pennsylvanian interglacials in comparison to pCO₂ during glacial periods (161 to 299 ppm ($-\mathrm{96}/+\mathrm{269}$ ppm)), with interglacial durations on the order of thousands to tens of thousands of years given the inferred eccentricity scale duration of the glacial–interglacial cycles (Horton et al., 2012; Montañez et al., 2016; Fielding et al., 2020).

Overall, the new pCO₂ record documents declining CO₂ through the final 13 Myr of the Pennsylvanian into the earliest Permian, including a 2.5 Myr interval (307 to 304.5 Ma) of minimum CO₂ values (less than 400 ppm to as low as 200 ppm) in the Kasimovian (Figs. 2b, S1b). Declining pCO₂ in the late Carboniferous coincides with rising atmospheric pO₂ (Glasspool et al., 2015; Krause et al., 2018; Lenton et al., 2018); thus, O₂∕CO₂ ratios in the interval of minimum Pennsylvanian CO₂ are nearly 2 times those of present-day values (∼515; gray line in Fig. 3a). A 10 Myr CO₂ nadir (∼180 to <400 ppm) characterizes the first two stages (Asselian and Sakmarian; 298.9 to 290.1 Ma) of the early Permian, overlaps with the peak occurrence of glacial deposits in the LPIA (gray boxes in Fig. 2b; Soreghan et al., 2019), and defines a second interval of anomalously high O₂∕CO₂ ratios (up to 970 ppm; Fig. 3a). A subsequent long-term rise (∼17 Myr) in pCO₂ to peak values up to ∼740 ppm $(-\mathrm{190}/+\mathrm{258}\mathrm{ppm})$ defines the remainder of the early Permian and coincides with multiple episodes of extensive and long-lived volcanism (Fig. 2b; Torsvik et al., 2008; Zhai et al., 2013; Sato et al., 2015; Shellnutt, 2018; Chen and Xu, 2019). This pCO₂ rise is also coincident with a decline in O₂∕CO₂ to below present-day values (Fig. 3a, b).

Short-term intervals of rising or falling CO₂ in the LOESS trend, within dating uncertainties, coincide with a brief but acute glaciation in the Kasimovian and with repeated deglaciations in south-central Gondwana in the early Permian (Griffis et al., 2018, 2019), as well as with restructuring of marine and terrestrial ecosystems (Fig. 3b–d). The statistical significance of these short-term rises and falls in CO₂ was evaluated by analyzing the raw Monte Carlo estimates (10 000 model runs per data point shown on the LOESS trend) generated by the aforementioned CO₂ models (Breecker, 2013; Franks et al., 2014), from which the bootstrapped CO₂ estimates for 11 increments of short-term rise or fall were subsequently determined (Fig. 4a). The analysis of the Monte Carlo CO₂ estimates within these short-term intervals of rising or falling CO₂ indicates that 72.5 % to 100 % of the data confirm a visually observed increasing or decreasing trend (Fig. 4).

4 Discussion

4.1 Declining CO₂ through the main phase of the LPIA

Atmospheric CO₂ concentrations in the final 13 Myr of the Carboniferous (the Pennsylvanian portion of our record) are generally higher than those of the earliest Permian (Fig. 2b) and overall decline through the later part of the Carboniferous. Higher pCO₂ in the latter half of the Pennsylvanian is compatible with the hypothesized waning of large early-to-middle Pennsylvanian glaciers in the late Pennsylvanian (cf. Fielding et al., 2008), including widespread terminal deglaciation in a major glacial depocenter in south-central Gondwana (Parana Basin, Brazil) toward the close of the Carboniferous (Griffis et al., 2018, 2019). Declining pCO₂ toward a nadir in the earliest Permian is also consistent with a renewed increase in the geographic distribution of glacial deposits in Gondwana beginning in the late Pennsylvanian and peaking in the earliest Permian (Fig. 2b; Soreghan et al., 2019).

A tectonically driven increase in CO₂ consumption via a strengthening of the silicate weathering (“climate stabilizing”) negative feedback (Walker et al., 1981; Berner and Caldeira, 1997) has been proposed as the driver of the Pennsylvanian decline in pCO₂ (Goddéris et al., 2017).The strength of the negative feedback varies with the degree of “weatherability” (i.e., the susceptibility to weathering), which, in turn, is predominantly controlled by the intensity of the hydrologic cycle (precipitation and surface runoff), with further influence by surface temperature and vascular plants (Dessert et al., 2001; Donnadieu et al., 2004; West, 2012; Maher and Chamberlain, 2014; Caves et al., 2016; Ibarra et al., 2016). Uplift of the Central Pangean Mountains (CPM) through the Pennsylvanian would have increased weatherability in the tropics by inducing orographic precipitation and creating steeper slopes (Goddéris et al., 2017), thus providing a greater supply of fresh mineral surfaces and enhanced surface runoff with longer fluid travel paths (cf. Maher and Chamberlain, 2014). Consequently, CPM-induced increased weatherability and CO₂ consumption would have enhanced the global efficiency of weathering and created a tighter coupling between CO₂ and climate at this time (cf. Maher and Chamberlain, 2014; Caves et al., 2016).

Figure 5 Carboniferous through early Permian modeled (GEOCLIM) steady-state atmospheric CO₂ and seawater ⁸⁷Sr∕⁸⁶Sr for different surface areas of mafic rocks available for silicate weathering. In the model, maximum geographic extent and altitude (5000 m) of the CPM is reached in the Moscovian (320 Ma), with altitude decreasing to 3000 m at 290 Ma and 2000 m at 270 Ma. (a) Simulated (color lines) and proxy pCO₂ estimates (purple crosses, this study) are shown. Horizontal error bars on the colored lines represent the temporal uncertainty for simulated pCO₂ estimates. CO₂ thresholds for continental ice sheet initiation (360–340 Ma = 1120 ppm; 340–300 Ma = 840 ppm; 300–260 Ma = 560 ppm from Lowry et al., 2014) decrease in response to equatorward drift of Gondwana, favoring an overall reduction in ice sheet size through time. The reference surface area of outcropping mafic rocks (GEOCLIM REG) maintains steady-state atmospheric CO₂ below the ice initiation threshold from 350 to ∼304 Ma. Steady-state atmospheric CO₂ for a 2-fold, 3-fold, or 4-fold increase in outcropping area of mafic rocks remains below the ice initiation threshold (560 ppm) up to ∼300 Ma, crossing over at progressively later times in the early Permian. Threshold crossover of steady-state CO₂ at ∼290 Ma for a 4-fold increase in mafic rock exposure coincides with the termination of the 10 Myr CO₂ nadir (gray vertical bar; both panels). (b) Seawater ⁸⁷Sr∕⁸⁶Sr modeled for the same set of varying surface areas of outcropping mafic rocks and ⁸⁷Sr∕⁸⁶Sr values of well-preserved biogenic calcites (gray X's) and conodont bioapatites (green and blue X's). Horizontal error bars on the colored lines represent the temporal uncertainty for modeled seawater ⁸⁷Sr∕⁸⁶Sr.

We model a Himalayan-type mountain range (an analog for the CPM) using GEOCLIM, with the model parameterized such that 30 % of the alkalinity generated by silicate weathering originates from weathering of mafic rocks. This is referred to as the “reference” continental silicate weathering mineral assemblage or GEOCLIM REG. That modeling indicates steady-state CO₂ concentrations (blue lines in Fig. 5a) that are well below the middle-to-late Carboniferous (340 to 300 Ma) threshold for initiation of continental ice sheets (840 ppm; Lowry et al., 2014). A hypothesized primary influence of the CPM on CO₂ consumption through increased weatherability is further supported by the coincidence of modeled seawater and marine proxy ⁸⁷Sr∕⁸⁶Sr values that define a plateau of peak radiogenic values that is sustained for 15 Myr of the late Carboniferous (318 to 303 Ma; Fig. 5b). The proxy-based seawater ⁸⁷Sr∕⁸⁶Sr plateau has been long interpreted to record exposure and weathering of uplifted and metamorphosed crustal rocks of the CPM that had radiogenic Sr isotope compositions (Chen et al., 2018, 2020). Additionally, the burial of substantial organic matter as peat in swamp environments prone to preservation (ultimately as coal) during the Pennsylvanian would have partitioned global CO₂ consumption between silicate weathering and organic carbon burial, further driving steady-state pCO₂ lower (D'Antonio et al., 2019; Ibarra et al., 2019). Our modeling, however, assumes a constant pre-Hercynian solid Earth degassing through the study interval, and it does not account for increased magmatic CO₂ during Hercynian arc-continent collision and potential widespread eruptive volcanism in the late Carboniferous (Soreghan et al., 2019), both of which could have increased steady-state CO₂.

Short-term fluctuations in pCO₂ are superimposed on the long-term decline through the latter portion of the Carboniferous. These short-term fluctuations have been confirmed as statistically significant (99.9 % to 100 % of estimates; Fig. 4b–d) and coincide with major environmental and biotic events. The brief interval of minimum pCO₂ (an average of ∼300 ppm but as low as 180 ppm) in the late Carboniferous (Kasimovian stage, 307 to 304.5 Ma; Fig. 3b) coincides with a short-lived but acute glaciation (306.5 to 305 Ma) recorded by prominent valley incision and large-scale regression recorded by cyclothemic successions in the US Appalachian Basin and Midcontinent basin, as well as Donets Basin, Ukraine (Belt et al., 2011; Eros et al., 2012; Montañez et al., 2016). Significant and repeated restructuring of wetland forests throughout tropical Euramerica, involving quantitative changes in floral composition and dominance, occurred during this 2.5 Myr pCO₂ minimum (and O₂∕CO₂ maximum; Fig. 3a–c). Before the short-term pCO₂ low, Euramerican tropical forests had expanded to their maximum aerial extent ($\ge \mathrm{2}\times {\mathrm{10}}^{\mathrm{6}}$ km²) under CO₂ concentrations of ∼500 ppm (Moscovian stage, Fig. 3b). The aerial extent of these forests dropped by half (green X in Fig. 3c; Cleal and Thomas, 2005) coincident with the decline in pCO₂ and near doubling of O₂∕CO₂ (Fig. 3a, b). Moreover, within this pCO₂ low (Fig. 3b), arborescent lycopsids of the wetland forests went extinct throughout Euramerica (white X in Fig. 3c) and seasonally dry tropical floras shifted from cordaitalean dominance to walchian dominance (∼307–306.8 Ma; Fig. 3c; DiMichele et al., 2009; Falcon-Lang et al., 2018). These restructuring events occurred at or proximal to CO₂ falling below 400 ppm, supporting a previously hypothesized but untested CO₂ threshold for the Pennsylvanian ecologic turnovers (Fig. 3b, c; Beerling et al., 1998; Beerling and Berner, 2000; Montañez et al., 2016). In the oceans, foraminiferal diversity decreased substantially during the Kasimovian pCO₂ low with the loss of ∼200 species (∼58 % of all taxa; first gray bar in Fig. 3d; Groves and Yue, 2009), presumably due to decreasing seawater temperatures.

The interval of CO₂ minima was terminated by a rapid rise across the Kasimovian–Gzhelian boundary (303.7 Ma) to CO₂ concentrations above 600 ppm (Figs. 2b; S1b). The short-term interval of elevated pCO₂ (304 to 302.5 Ma) is coincident with a ∼1.5 ‰ decline in seawater δ¹³C (Grossman et al., 2008; Chen et al., 2020), which is compatible with a decline in the CO₂ sink provided by terrestrial organic C (peats) burial (gray bar in Fig. 2b) and∕or a peak in pyroclastic volcanism between ∼310 and 301 Ma (Soreghan et al., 2019). This period of increased pCO₂ overlaps with the Alykaevo Climatic Optimum (orange bar in Fig. 3c), defined by the invasion of tropical Euramerican vegetation into the Rufloria-dominated midlatitude Angaran floral province (Cleal and Thomas, 2005). Terminal deglaciation in south-central Gondwana (Parana Basin, Brazil), U–Pb dated to between ∼302 and 298 Ma (Cagliari et al., 2016; Griffis et al., 2018), may have been linked to the late Pennsylvanian interval of elevated CO₂, although this requires further testing (Figs. 2b, 3b). Conversely to the Kasimovian CO₂ low, a significant change in global diversity of foraminifera involving a doubling of species occurred during this subsequent period of elevated CO₂ and presumed increase in seawater temperatures (black bar in Fig. 3d; Groves and Yue, 2009).

4.2 An early Permian CO₂ nadir

Atmospheric pCO₂ dropped substantially across the Carboniferous–Permian boundary (i.e., 298.9 Ma) to a 10 Myr interval (300–290 Ma) of the lowest concentrations (175 to <400 ppm) of the 40 Myr record (Fig. 2b). The CO₂ nadir, which spans the Asselian and Sakmarian stages, coincides with renewed glaciation and maximum ice sheet extent, marking the apex of LPIA glaciation (Fig. 2b; Fielding et al., 2008; Isbell et al., 2012; Montañez and Poulsen, 2013; Soreghan et al., 2019), as well as with a large-magnitude eustatic fall archived in paleotropical successions worldwide (Koch and Frank, 2011; Eros et al., 2012). Widespread glacial expansion temporally linked to this interval of lowest overall pCO₂ argues for CO₂ as the primary driver of glaciation rather than recently proposed mechanisms, such as the influence of the closing of the Precaspian Isthmus (Davydov, 2018) or a decrease in the radiative forcing resulting from increased atmospheric aerosols by explosive volcanism at this time (Soreghan et al., 2019). The very low greenhouse radiative forcing associated with this low CO₂ interval would have been amplified by 2.5 % lower solar luminosity (Crowley and Baum, 1992), reduced transmission of shortwave radiation (Poulsen et al., 2015) by the high pO₂ atmosphere of the early Permian (Krause et al., 2018; Lenton et al., 2018), and by increased atmospheric aerosols at this time (Soreghan et al., 2019).

The 10 Myr pCO₂ nadir raises a paradox as to what was the primary CO₂ sink(s) at the time given that the CO₂ sinks of the Pennsylvanian were no longer prevalent. This paradox reflects the waning denudation rates of the CPM by the early Permian (Goddéris et al., 2017), intensifying pantropical aridification, possibly driven by increasing continentality (yellow-to-red bar in Fig. 3c; DiMichele et al., 2009; Tabor et al., 2013) and the demise of the wetland tropical forests and associated loss of peats before the close of the Carboniferous (black-to-gray bar in Fig. 2b; Hibbett et al., 2016). In turn, surface runoff would have been inhibited and the supply of fresh silicate minerals decreased, thus lowering overall weatherability. Atmospheric CO₂ under the influence of these environmental factors should have equilibrated in the earliest Permian at a new higher steady-state level, even if solid Earth degassing did not increase (cf. Gibbs et al., 1999), thus raising a paradox. If volcanism was increasing by this time (Fig. 2b and associated references; Soreghan et al., 2019), then this paradox is even greater.

Two mechanisms have the potential to resolve this paradox. The first, referred to as a delayed climate-controlled capacitor (Joshi et al., 2019), leads to a multi-million-year delay between the timing of peak orogenic uplift and maximum chemical weathering potential and CO₂ drawdown due to substantial differences in chemical weathering rates during the different phases of an orogenic cycle. In their study, the highest intensity of chemical weathering and capacity for CO₂ consumption occurs when mountains have been somewhat denuded rather than during peak uplift, reflecting the disproportionate influence of runoff temperature over hydrology and erosion on weathering potential. Notably, the coupled climate and geochemical modeling of the late Paleozoic ice age by Joshi et al. (2019) yields an evolution of simulated pCO₂ over the period of uplift and denudation of the CPM that corresponds to both absolute CO₂ concentrations and the magnitude of change over this period (∼320 to 290 Ma). That said, in Joshi et al. (2019), the physical erosion parameter is not dependent on climate but rather is defined by the prescribed altitude. Thus, physical erosion is an external forcing in their model. The absence of runoff dependence for physical erosion (as is the case for GEOCLIM) and the strong dependence of weathering on temperature may be the trigger for their simulated delay between maximum uplift and the highest intensity of CO₂ consumption by silicate weathering. In GEOCLIM, the dependence of physical erosion on runoff does not allow for millions of years delay between maximum uplift of the CPM and lowest simulated pCO_2. Further study is needed to interrogate the influence of this approach on the results.

The second mechanism, proposed here, is a substantial shift in the ratio of mafic-to-granite rocks available for weathering from the latest Carboniferous to the early Permian. This reflects the doubling or greater increase in weatherability of mafic mineral assemblages over granitic assemblages (Gaillardet et al., 1999; Dessert et al., 2003; Ibarra et al., 2016), thus enhancing weathering efficiency and CO₂ drawdown and creating a tighter coupling between CO₂ and climate. In turn, with tighter coupling between CO₂ and climate, the global silicate weathering flux needed to maintain homeostatic balance in the carbon cycle for a given scenario can be attained at a lower pCO₂ level.

Macdonald et al. (2019) hypothesized that increased weatherability provided by the exhumation of ophiolites along the ∼10 000 km long Hercynian arc-continent suture zone, primarily situated in the paleotropics, was capable of lowering pCO₂ below the ice initiation threshold in the Carboniferous (i.e., Pennsylvanian), thus instigating the late Paleozoic ice age. Here, we used the GEOCLIM model to, first, interrogate this Carboniferous hypothesis further and, second, to evaluate the potential of increased weatherability provided by increasing the ratio of outcropping mafic rocks to granite rocks available for weathering, as the predominant driver of the early Permian CO₂ nadir. Figure 5 illustrates the influence of a successive increase in the surface area of outcropping mafic rocks beginning with the reference continental silicate mineral assemblage (GEOCLIM REG), which was used to evaluate the influence of Pennsylvanian uplift of the CPM to an up to 4-fold increase in the outcropping of mafic rocks. In the GEOCLIM context, the weathering of mafic rocks is dependent on the surface of each grid cell (and on the associated local runoff and air temperature) multiplied by a calibration constant. Increasing the exposure area of mafic rocks is mathematically equivalent to multiplying by the calibration constant.

Between 300 and 290 Ma, when predominant Pennsylvanian CO₂ sinks were lost (terrestrial organic C burial) or waning (decreased precipitation and denudation rates of the CPM), modeled steady-state atmospheric CO₂ was maintained at or below the CO₂ threshold for initiation of continental ice sheets (560 ppm; Lowry et al., 2014) when the surface area of outcropping mafic rocks was greater than 2-fold that of GEOCLIM REG (Fig. 5a). Conversely, steady-state CO₂ rises well above the glacial threshold (to 3500 pm) for the reference continental silicate rock assemblage (Fig. 5a). Although volcanism remained geographically extensive through the 10 Myr CO₂ nadir (Soreghan et al., 2019), the impact on atmospheric CO₂ would have been short lived (≤10⁵ kyr; Lee and Dee, 2019) and eclipsed on the longer term by the increased weatherability provided by increased exposure of mafic rocks along the Hercynian arc-continent suture zone, lowering steady-state CO₂ to, potentially, pre-volcanism levels (cf. Dessert et al., 2001).

Independent evidence for a substantial shift in the partitioning of silicate weathering to more mafic mineral assemblages in the earliest Permian exists in the late Paleozoic proxy-based seawater Sr isotope record, which documents a rapid (0.000043 Myr⁻¹) and near-linear decrease in seawater ⁸⁷Sr∕⁸⁶Sr beginning in the latest Carboniferous (∼303 Ma) and continuing into the middle Permian (Fig. 5b; Chen et al., 2020). The simulated trends in seawater ⁸⁷Sr∕⁸⁶Sr for GEOCLIM REG (blue line on Fig. 5b) through a 2- to 4-fold increase in the area of exposed mafic rocks capture the rapid rise through the upper Carboniferous to peak values in the latter half of the Pennsylvanian and subsequent decline through the early Permian. The rapid rate of decline in proxy ⁸⁷Sr∕⁸⁶Sr values post-300 Ma, however, is best bracketed by simulated ⁸⁷Sr∕⁸⁶Sr for a 2- to 4-fold increase in mafic rock exposure. Moreover, the best fit of the simulated trends to the geochronologically well-constrained bioapatite data (blue and green crosses in Fig. 5b) suggests a progressive increase in mafic-to-granite ratio through the 10 Myr CO₂ nadir. This finding, together with the hypothesized need (the aforementioned second mechanism) for minimally a 4-fold increase in mafic-rock outcropping in order to maintain CO₂ concentrations below the ice initiation threshold for a sustained period longer than that of hysteresis (i.e., throughout the interval of minimum CO₂ and apex of glaciation; Fig. 5), supports a change in weatherability through time. That is, the modeling results argue for a substantial increase in weatherability from the Carboniferous to early Permian driven by a compositional shift in outcropping rocks available for weathering from that of the reference continental silicate assemblage to one of higher mafic-to-granite ratio.

If peak ophiolite exhumation and maximum CO₂ consumption by their weathering occurred in the late Carboniferous, thus initiating the LPIA (∼330 to 300 Ma) as has been suggested (Table S1 of Macdonald et al., 2019), then our modeling results suggest that a substantial increase in solid Earth degassing rate at this time would have been necessary. In our simulation, increasing the surface area of outcropping mafic rocks (2- to 4-fold) during the Pennsylvanian results in steady-state atmospheric CO₂ levels approaching “Snowball Earth” conditions given other operating influences on weatherability and CO₂ sequestration at the time and no change in degassing rate (Fig. S6). Such conditions are not compatible with proxy-inferred moderate surface conditions of the late Carboniferous (Montañez and Poulsen, 2013) and the radiation of forest ecosystems throughout the tropics (DiMichele, 2014). Rather, we hypothesize that the sustained CO₂ nadir and expansion of ice sheets in the first 10 Myr of the Permian record a major reorganization of the predominant factors influencing weatherability in the tropics across the Carboniferous–Permian transition, in particular a substantial shift in the ratio of mafic-to-granitic rocks available for weathering.

Similar to the short-term fluctuations superimposed on the later Carboniferous long-term decline in CO₂, two statistically significant (94 % to 100 % in Fig. 4e–h) short-term increases in pCO₂ are superimposed on the early Permian nadir (Fig. 3b). The first (298 to 296 Ma) coincides, within age uncertainty, with a major deglaciation event in the Karoo (southern Africa) and Kalahari (Namibia) basins of south-central Gondwana (296.41 Ma $+\mathrm{0.27}/-\mathrm{0.35}$ Ma; Griffis et al., 2019). The second short-term rise in pCO₂ (294.5 to 292.5 Ma) overlaps with the onset of widespread ice loss in several southern Gondwanan ice centers (Fig. 2b; Soreghan et al., 2019). This CO₂–deglaciation link suggests that continental ice stability in the early Permian dropped substantially when pCO₂ rose above ∼300 to 400 ppm and thus raises the question as to whether the ice sheet CO₂ threshold was even lower than modeled (560 ppm; Lowry et al., 2014) during the earliest Permian.

Notably, the early Permian pCO₂ minimum (Figs. 2b, S1b), as well as the associated anomalously high O₂∕CO₂ ratios (700 to 960; Fig. 3a), is an interval of major ecosystem changes. The geologically rapid and large-magnitude drop in pCO₂ prior to and across the Carboniferous–Permian boundary interval (Figs. 2b, S1b) coincides with a permanent shift in plant dominance from swamp-community floras to seasonally dry vegetation (black X in Fig. 3c). That shift in plant dominance has been long attributed to intensification of an aridification trend that began in the mid-Pennsylvanian (yellow-to-red bar in Fig. 3c; DiMichele et al., 2009; Tabor et al., 2013). Analogous to the vegetation turnover and extinction during the Pennsylvanian CO₂ minimum, this permanent shift to seasonally dry vegetation is coincident with the earliest Permian drop in pCO₂ to concentrations below 400 ppm, suggesting a possible ecophysiological advantage of these plants over the wetland floral dominants that they replaced (Fig. 3a–c; cf., Wilson et al., 2017). The high water use efficiency of the seasonally dry plants would have made them water-stress tolerant and, in turn, would have amplified aridification through a modeled ∼50 % decrease in canopy-scale transpiration (Wilson et al., 2017, 2020). The extreme habitat restriction of wetland floras was particularly consequential for tetrapods, leading to the acquisition of terrestrial adaptions in crown tetrapods and the radiation and eventual dominance of dryland-adapted amniotes, possibly shaping the phylogeny of modern terrestrial vertebrates (Fig. 3c; Pardo et al., 2019).

Moreover, the CO₂ decline at the Carboniferous–Permian boundary into the 10 Myr nadir and associated peak in O₂∕CO₂ also corresponds to the evolution and radiation of glossopterids and gigantopterids (McLoughlin, 2011; Zhou et al., 2017), with increasing vein density in the former (Fig. 3a–c; Srivastava, 1991). These plant groups had complex angiosperm-like venation (Melville, 1983; Srivastava, 1991), with gigantopterids having the only known pre-Cretaceous vessels in their stems (Li et al., 1996). The increased hydraulic capacity provided by these morphological characteristics would have conferred a significant ecological advantage to these plants under the low CO₂, high O₂, and elevated aridity conditions in which they evolved (cf. Gerhart and Ward, 2010; de Boer et al., 2016). In the oceans, a marked collapse in foraminiferal diversity with a notable fall in species to a minimum from a Pennsylvanian zenith (425 to 110 species; Fig. 3d, e; Groves and Yue, 2009) spanned the 10 Myr pCO₂ nadir, analogous to the diversity drop during the Pennsylvanian low-CO₂ interval.

4.3 CO₂-forced demise of the LPIA

The 10 Myr CO₂ nadir terminated at 290 Ma with the onset of a protracted CO₂ rise that persisted to the highest levels of the record (∼740 ppm ($-\mathrm{190}/+\mathrm{258}$)) by the close of the early Permian (Fig. 2b). The onset of this protracted CO₂ rise overlaps with initiation of a period of large-magnitude magmatism (red bars in Fig. 2b). Widespread volcanism began around 297.4 Ma (±3.8 Ma) in northern Europe (Skagerrak-centered large igneous province), extending well into Germany (Rotliegend) (Torsvik et al., 2008; Käßner et al., 2019). The multistage Tarim magmatic episodes in China (292–272 Ma; with peaks at ∼290 and 280 Ma; Fig. 2b; Chen and Xu, 2019), were likely associated with large-magnitude CO₂ emissions given that the magma, which distributed basalt (400 m thick) over a 2.5×10⁵ km² region (Yang et al., 2013), intruded a thick succession of early Paleozoic marine carbonates (Gao et al., 2017). The Panjal Traps, NW India (289±3 Ma; Shellnutt, 2018), and the compositionally similar Qiangtang dikes (283±2 Ma; Fig. 2b; Zhai et al., 2013), albeit relatively small in extent, were an additional potential volcanic CO₂ source, along with contemporaneous volcanism in Oman. Furthermore, protracted Choiyoi volcanism, which began at 286.5±2.3 Ma (Sato et al., 2015) and continued over ∼39 Myr in western Argentina, may have contributed substantial pulses of greenhouse gases in the early Permian (Spalletti and Limarino, 2017). Once each magmatic episode waned, however, the mafic-dominated magmatic deposits would have served as longer-term regional sinks, leading to increased global CO₂ consumption (cf. Lee et al., 2015). Thus, for steady-state CO₂ to have increased through the remainder of the early Permian, the relative influence of CO₂ inputs must have outpaced that of these, and other, outputs (CO₂ sinks).

Our modeled (GEOCLIM) steady-state CO₂ for a 4-fold increase in outcropping of mafic rocks surpasses the ice-sheet initiation threshold at the termination of the CO₂ nadir (∼290 Ma; red line in Fig. 5a) despite no change in solid Earth degassing. The fact that low CO₂ concentrations could no longer be maintained, despite a 4-fold increase in mafic rock exposure, reflects overall intensifying aridification, denudation of the CPM, and a shift from dense forests to savanna-like vegetation in Euramerica at this time. However, given that the magmatic CO₂ flux likely increased already by the earliest Permian (summarized in Fig. 2b), our model results indicate that maintaining low steady-state CO₂ concentrations during the earliest Permian 10 Myr CO₂ nadir would have required an increasingly greater proportion of mafic rock weathering over the reference continental silicate mineral assemblage of the Pennsylvanian, possibly well beyond a 4-fold increase.

A CO₂-forced demise of the late Paleozoic ice age after 290 Ma is supported by the loss of continental ice from the main ice depocenters in south-central Gondwana by $\mathrm{282.17}+\mathrm{0.32}/-\mathrm{0.44}$ Ma (Griffis et al., 2018, 2019) and a 6-fold drop in documented glacial deposits overall between the Sakmarian and Artinskian stages (Fig. 2b; Soreghan et al., 2019). The long-term CO₂ rise through the remainder of the early Permian coincided with substantial marine and terrestrial ecosystem perturbation (Fig. 3b–d; Chen and Xu, 2019). In the marine biosphere, the uniformly low rates of global macroevolution in marine organisms (brown bar in Fig. 3d) were reversed, and narrowly adapted and distributed genera reappeared, thus restoring marine ecosystems to their pre-LPIA rates (Stanley and Powell, 2003; Powell, 2005). Pennsylvanian rugose corals (pink bar in Fig. 3d) underwent a major turnover in composition to those that dominated until the end-Permian extinction, and cold-adapted marine bivalves and brachiopods turned over to warm-adapted forms synchronous with the onset of the long-term increase in pCO₂ (290.1 Ma; blue to red bar across the Sakmarian–Artinskian boundary in Fig. 3d; Wang et al., 2006; Clapham and James, 2008). On land, the loss of pelycosaur families (three in the late Artinskian and four in the early Kungurian; Kemp, 2006) coincided with CO₂ sustained at >500 ppm. By the close of the Kungurian and the time of highest CO₂ (740 ppm), basal synapsids largely disappeared and were replaced by more derived therapsids, tetrapod diversity decreased significantly (Benton, 2012; McGhee, 2018), plant extinction rates reached a level comparable to that associated with the extinction of arborescent lycopsids in the early Kasimovian (Cascales-Miñana et al., 2016), and extinction and origination rates increased in fishes (Friedman and Sallan, 2012).

5 Conclusions

Glacial–interglacial climate cycles and large-scale glacioeustasy, as well as repeated ecosystem change analogous to that of the Pleistocene, characterized Earth's penultimate icehouse period in the late Paleozoic. The dynamic glaciation history of this icehouse period (the late Paleozoic ice age – LPIA) came to a close by the end of the early Permian with turnover to permanent greenhouse conditions. Thus, improved constraints on how atmospheric pCO₂ evolved during the LPIA and its subsequent demise is crucial for better understanding the role of greenhouse gas forcing on Earth system processes during this time. The new and age-recalibrated pCO₂ reconstruction presented here for a 40 Myr interval (∼313 to 273 Ma) of the late Paleozoic substantially refines existing Permian CO₂ estimates and provides perhaps the highest temporal resolution protracted pCO₂ record prior to the Cenozoic. The multiproxy record confirms the previously hypothesized CO₂–glaciation linkage, including documenting the coincidence of a 10 Myr period of minimum pCO₂ with inferred maximum ice extent during the earliest Permian. A long-term decline in pCO₂ through the late Carboniferous period of glaciation, culminating in the earliest Permian CO₂ nadir, lends support for a previously modeled progressive decrease in the CO₂ threshold for continental ice sheets through the LPIA.

Our new pCO₂ record provides the first stomatal-based evidence for elevated (up to 700 ppm) atmospheric CO₂ concentrations during short-term (10⁴-year) interglacials. Together with new O₂∕CO₂ estimates of similar temporal resolution to pCO₂, the new atmospheric trends indicate a close temporal relationship to repeated ecosystem restructuring in the terrestrial and marine realms. In terrestrial ecosystems, the appearance and∕or rise to dominance of plants with physiological and anatomical mechanisms for coping with CO₂ starvation and marked aridity correspond to drops in CO₂ below 400 ppm (as low as ∼180 ppm), and O₂∕CO₂ ratios nearly doubled those of late Paleozoic background values. Similarly, decreasing rates of macroevolution and diversity in the low-latitude oceans corresponded to falling CO₂ to below 400 ppm. These CO₂–ecosystem relationships lead us to hypothesize that 400 ppm was an important threshold for ecosystem resilience during the late Paleozoic.

Modeling of steady-state pCO₂ during the late Paleozoic using an intermediate-complexity climate–carbon cycle model (GEOCLIM) and comparison to the new multiproxy CO₂ record provides new insight into the relative influences of uplift of the Central Pangean Mountains, intensifying aridification, and increasing mafic to granite ratio of outcropping rocks on the global efficiency of CO₂ consumption and secular change in steady-state pCO₂ through the late Paleozoic. The simulations confirm that, for the Carboniferous, enhanced weatherability and CO₂ consumption provided by the influence of uplift on surface hydrology and erosion could have lowered atmospheric pCO₂ well below the threshold for ice sheet initiation. Increasing the availability of mafic rocks for weathering drives CO₂ levels toward Snowball Earth conditions in the Carboniferous. Conversely, a substantial increase (up to a 4-fold increase) in the surface outcropping of mafic rocks over those modeled for the Carboniferous is needed to maintain the 10 Myr CO₂ nadir in the earliest Permian and is compatible with maximum exhumation of the Hercynian orogenic belt at this time, as well as with a rapid decline in proxy-based seawater ⁸⁷Sr∕⁸⁶Sr. Although these findings support the hypothesis of atmospheric pCO₂ response to uplift of the CPM as the primary driver for Carboniferous initiation of the LPIA (Goddéris et al., 2017), they argue for a major reorganization of the predominant surface factors influencing weatherability in the tropics across the Carboniferous–Permian in order for pCO₂ values to have fallen to well below 300 ppm. The demise of the LPIA was greenhouse gas forced, reflecting the increasing importance of magmatic degassing and likely decreased weathering efficiency driven by intensifying aridification, denudation of the CPM, and the loss of the wetland forests throughout tropical Euramerica.