Marine anoxia linked to abrupt global warming during Earth's penultimate icehouse.

Piecing together the history of carbon (C) perturbation events throughout Earth’s history has provided key insights into how the Earth system responds to abrupt warming. Previous studies, however, focused on short-term warming events that were superimposed on longer-term greenhouse climate states. Here, we present an integrated proxy (C and uranium [U] isotopes and paleo CO2) and multicomponent modeling approach to investigate an abrupt C perturbation and global warming event (∼304 Ma) that occurred during a paleo-glacial state. We report pronounced negative C and U isotopic excursions coincident with a doubling of atmospheric CO2 partial pressure and a biodiversity nadir. The isotopic excursions can be linked to an injection of ∼9,000 Gt of organic matter–derived C over ∼300 kyr and to near 20% of areal extent of seafloor anoxia. Earth system modeling indicates that widespread anoxic conditions can be linked to enhanced thermocline stratification and increased nutrient fluxes during this global warming within an icehouse.

Observations and climate models indicate that the dissolved oxygen inventory of the modern ocean is decreasing, with temperature-driven decline in oxygen solubility being a key driver (1). A decline in ocean dissolved oxygen, expressed as an expanded oxygen minimum zone (OMZ), will negatively impact marine ecosystems, leading to significant loss of biodiversity in the ocean (2). This is of significant concern for the world’s largest fisheries situated in the most productive areas of global oceans, as these regions are particularly susceptible to ocean deoxygenation (3, 4). Substantial uncertainty in estimating the extent of deoxygenation over the upcoming millennia, however, drives the current focus on understanding past episodes of ocean deoxygenation.

Empirical constraints on the magnitude of ocean deoxygenation during climate perturbations come predominantly from the Quaternary glacial–interglacial transitions (5, 6) or early Cenozoic rapid warming events, in particular the Paleocene–Eocene Thermal Maximum (PETM) event (7, 8). The temporal scales of warming and deoxygenation of these events differ by an order of magnitude (104 vs. 105 y). Constraints on ocean circulation and biogeochemical cycles across warming events in the Quaternary are more robust than in Earth’s deep past (5, 9). On the other hand, Quaternary partial pressure of CO2 (pCO2) and temperature shifts were gradual and the overall perturbations small in magnitude (10) relative to predicted changes for the next millennia or two. Although changes in pCO2 during the early Cenozoic warming events (foremost, the PETM) were larger in magnitude and more rapid than carbon (C) perturbations of Quaternary glacial–interglacial transitions (11, 12), they occurred under a background greenhouse climate state characterized by high baseline atmospheric pCO2 (∼1,000 ppm). Other periods of pre-Cenozoic C perturbations (13, 14), such as the Cretaceous and Jurassic (Toarcian) ocean anoxic events (OAEs) (15, 16), and the end-Triassic (17) and the end-Permian mass extinction events (18), also occurred during background greenhouse climates (19, 20). These greenhouse OAEs have provided constraints on and insights into how to model climate change and marine redox evolution. To date, the degree of deoxygenation and spread of anoxic conditions with warming in a glacial state is relatively unexplored.

Here, we provide a perspective on global warming–induced ocean deoxygenation by documenting a 105-y C-perturbation event associated with widespread oceanic anoxia, superimposed on the late Paleozoic glacial climate state. This past icehouse (with main episode between ∼340 and 290 Ma) existed under atmospheric CO2 levels comparable to that of the past few million years (21) and was a period of dynamic and widespread glaciation in the Southern Hemisphere (Gondwana). As such, it provides unique constraints on deoxygenation with warming under a background glacial state, albeit with different paleogeographic boundary and marine and ecosystem conditions compared with those of the younger greenhouse periods of C perturbation.

Results and Discussion

Global Warming and C Emission across the Kasimovian–Gzhelian Boundary.

Several lines of evidence confirm a previously hypothesized long-term global warming in the late Pennsylvanian (21–23) and reveal an abrupt warming (105-y scale) across the Kasimovian–Gzhelian boundary (KGB; ∼304 Ma; Fig. 1). Foremost, a multiproxy-based reconstruction of atmospheric pCO2 demonstrated a pronounced increase in pCO2 from a baseline of ∼350 ppm to ∼700 ppm over ∼300 kyr just prior to the KGB (21, 24). This abrupt doubling of pCO2 followed on the heels of an ∼1.5-Myr nadir in pCO2 and acute glaciation. Second, sea-surface temperatures (SSTs) inferred from well-preserved brachiopod calcite oxygen isotopes (25) increased from ∼25 °C to ∼29 °C across this boundary (Fig. 1). Third, there was a widespread rise in sea level over several million years of the Middle to Late Pennsylvanian, culminating in maximum transgression across the KGB (26), consistent with a major loss in continental ice volume. Finally, an inferred warming, the Alykaevo Climatic Optimum, saw the migration of paleo-tropical flora into the Northern Hemisphere midlatitudes (22). Notably, the magnitude of change in pCO2 and carbonate δ13C (δ13Ccarb) across the KGB warming (Fig. 1) is substantially larger than the reconstructed glacial–interglacial variability (21), indicating that the KGB was an unusually pronounced interglacial.

Fig. 1.

Geochemical proxies and biodiversity records implying a significant perturbation to the C cycle and warming, and a nadir in biodiversity across the KGB. (A) Atmospheric pCO2 reconstruction based on multiple proxies with a locally estimated scatterplot smoothing (LOESS) trendline (dark gray; 15% smoothing) and bootstrapped errors (shaded; 75% CI) from ref. 24. (B) δ13C of calcitic brachiopods (squares) from ref. 25 and marine bulk carbonates (circles) from refs. 27 and 28 with a LOESS trendline (dark gray; 20% smoothing) and bootstrapped errors (shaded; 2 SD). (C) Oxygen isotopes (δ18O) of calcitic brachiopods (gray) from ref. 25 with a LOWESS trendline (dark gray; 30% smoothing) and bootstrapped errors (shaded; 2 SD), and calculated SSTs (blue), based on ref. 79. (D) Biodiversity reconstruction based on global compilation of foraminifera (black dashed line) from ref. 56 and all marine species of China from the Geobiodiversity Database (blue line) (57). Mid. Penn., Middle Pennsylvanian; P, Permian; A, Asselian; ppmv, parts per million volume; V-PDB, Vienna Pee Dee Belemnite.

Here, we present paired δ13Ccarb and organic carbon isotope (δ13Corg) records that span 4 Myr of the latest Carboniferous from South China (Fig. 2), providing insights into C cycling during the KGB warming event. Our records were developed from two open-water slope carbonate successions—the Naqing and Narao sections ()—both which were the targets of previous high-resolution δ13Ccarb studies (27, 28). The Naqing (N 25° 14′ 40″, E 106° 29′ 26″) and Narao (N 25° 24′ 39″, E 106° 36′ 25″) successions (∼22 km apart) were deposited on carbonate slopes of the Qian-Gui Basin, now located in the Luodian region, Guizhou Province, in southern China (). A compilation of bulk δ13Ccarb from the two South China successions captures a preboundary, negative C isotope excursion (∼2‰, from values of 4 to 2‰; Fig. 2), which has been documented worldwide, albeit with relatively low sampling resolution and low-resolution age constraints for most successions (25, 29). The negative isotopic excursion is also reflected in the δ13Corg records of bulk sediment from both South China successions (∼4‰, from values of −20 to −24‰; Fig. 2). Both δ13Ccarb and δ13Corg returned to background values within ∼300 kyr following the KGB. The difference in magnitude of the negative excursion between δ13Ccarb and δ13Corg is consistent with a rise in atmospheric pCO2 (30). Collectively, the coupled negative excursions of δ13Ccarb and δ13Corg in the Naqing and Narao successions indicate a perturbation to the C cycle across the KGB event.

Fig. 2.

Paired δ238U and δ13C data from carbonate slope successions (Naqing and Narao) in South China exhibiting pronounced negative excursions in δ238U and δ13C immediately below the KGB. (A and D) Chrono- and biostratigraphy and sedimentologic logging. (B and E) Carbonate δ238U data (±2 SD) with a locally estimated scatterplot smoothing (LOESS) trendline (gray; 20% smoothing). (C and F) Bulk δ13Ccarb (blue) from refs. 27 and 28 with a LOESS trendline (gray; 10% smoothing) and bulk rock organic matter δ13C (orange). For all panels, the bootstrapped errors are indicated by the shaded region (95% CI). I., Idiognathodus.

We used the global C-cycle model Long-Term Ocean-Atmosphere-Sediment Carbon Cycle Reservoir (LOSCAR) (31) and paleosol- and stomatal-based pCO2 estimates (21, 24) across the KGB interval to quantify the amount of C emitted, the duration of C injection, and the δ13C of C released that contributed to the doubling of pCO2 and the negative C-isotope excursion. A similar methodology has been used to explore the source of the C flux during the PETM (31, 32), and estimates from the LOSCAR model generally match models of greater complexity (12). To match the baseline δ13Ccarb (∼3.8‰) and pCO2 (∼350 ppm) before the KGB, we modified the default modern boundary settings (icehouse) of the LOSCAR model by adjusting the steady-state pCO2 level, the initial δ13C value of weathered carbonate, and the ratio of shelf carbonate precipitated relative to deep carbonate rain (Materials and Methods).

To quantify the amount, duration, and isotopic composition of the C emitted, we performed Monte Carlo resampling for these parameters and ran 432 LOSCAR models through the KGB (each model with a unique parameter combination). Modeled δ13C and pCO2 values were compared with the observed δ13Ccarb time series and empirical estimates of pCO2 to constrain the range of C-emission trajectories (Materials and Methods and ). Based on data-model comparisons, we estimate the total C released to range from 5,000 Gt (or 5,000 Pg) to 10,000 Gt, with 9,000 Gt as the most likely estimate (Fig. 3 and ). These results also point toward a δ13C signature of the emissions of approximately −19‰ released over ∼300 kyr (for sensitivity test, see ).

Fig. 3.

Global C-cycle (LOSCAR) model results (best scenario runs) and U mass-balance modeling results across the KGB. (A) LOSCAR modeling (blue line) of atmospheric pCO2 data (locally estimated scatterplot smoothing (15% smoothing; gray line). (B) LOSCAR modeling (blue line) of average surface seawater δ13C inferred from the Naqing succession, given its tighter age constraints (). (C) Original carbonate δ238U data (gray symbols) from the Naqing succession (with a mean before and after the KGB for regression; gray line) and diagenetically corrected δ238U data (subtracting 0.27‰; blue line) used for U-cycle modeling. (D) U-model estimates of fanox, documenting a rapid increase in the areal extend of seafloor anoxia to near 20% immediately below the KGB. Color bar indicates frequency of the result from the Monte Carlo simulations. ppmv, parts per million volume.

The markedly negative isotopic value of the emitted C requires an organic source. One potential source of the isotopically light C is thermogenic methane released during initial intrusion of volcanics into organic-rich sediments (e.g., Carboniferous coals) associated with the Skagerrak-Centered (Jutland-Skagerrak) Large Igneous Provinces (33, 34). Better age constraints on these volcanic intrusions are necessary to robustly test this hypothesis. An alternative possible CO2 source is organic C released from tundra or permafrost peripheral to retreating glaciers driven by warming (21, 35, 36). Uncertainties exist in comparing the model output with the empirical time series, due to the limited pCO2 data for certain time intervals and limited age constraints, thereby preventing us from pinpointing the exact C source. Nonetheless, our results suggest that the KGB warming was linked to a massive release of organic matter–derived CO2.

Marine Anoxia during the KGB Warming.

To explore the extent of ocean deoxygenation that occurred during the KGB warming, we generated carbonate uranium (U) isotope data from the two South China successions (Fig. 2). U isotopic values (238U/235U, denoted as δ238U) of seawater recorded in marine carbonates are now widely used to reconstruct redox conditions of the global ocean (16, 37–42). Under anoxic or euxinic conditions, dissolved U(VI) is reduced to U(IV), during which the dominant 238U isotope is removed from ocean water to anoxic sediments, leaving seawater enriched in 235U. It is worth noting that U isotopes are likely a response proximal to the iron reduction zone of the redox ladder (43). Thus, U sequestration and fractionation are likely minimal in the upper portion of the anoxic window, potentially leading to variability in contemporaneous δ238U proxy records (44). The long residence time of U in the ocean (320 to 560 kyr for the modern era) suggests U is well mixed globally (45). Given that the U flux to the ocean is likely to have had a constant isotopic value of −0.3‰ ± 0.04‰ (2 SE) (46), the δ238U value of seawater records the relative proportions of U removed in oxic and anoxic settings globally—with negative δ238U excursions marking periods of expanded marine anoxia (39, 47).

We find that the carbonate δ238U records from the South China successions exhibit relatively high δ238U values of between −0.3‰ and 0‰ below and above the KGB interval and a distinct negative excursion across the KGB with minimum values of −0.54‰ in the Naqing and −0.43‰ in the Narao records (Fig. 2). Recent studies indicate that the δ238U values of shallow-water carbonates are consistently, by an average of ∼0.24 to 0.27‰, more positive than the contemporaneous seawater (47, 48), whereas pelagic carbonate δ238U values approximate the seawater signature, suggesting a diagenetic gradient may exist among depositional sites of different paleobathymetry (16, 47, 49, 50). The isotopic offset results from both fractionation of U isotopes during early diagenetic carbonate dissolution and recrystallization, and authigenic enrichment of isotopically heavy U(IV) in carbonates (47, 48). To explore this, we used both petrographic and geochemical criteria (e.g., Mn/Sr, Mg/Ca, and total organic C) to evaluate the potential for diagenetic alteration of our samples (). We found no petrographic or geochemical evidence of meteoric early diagenetic alteration or of late diagenetic recrystallization, calcite cements, and dolomitization ().

The variations between the Naqing and Narao records across the KGB can be accounted for by differences in local depositional and subseawater-sediment conditions (e.g., pore-water redox chemistry and sedimentation rate) between the sections. It is feasible that the less negative δ238U values across the KGB in the Narao section are due to the shallower water environment and greater authigenic enrichment of isotopically heavy U(IV). Similarly, the seawater δ238U values immediately below the KGB in the two sections might be shifted positively, given the relatively enriched total organic C contents in the bulk sediments (). The pattern of a relatively rapid decrease in δ238U to a nadir observed in the Naqing succession is expanded in, and is likely driven by higher sedimentation rates of, the Narao succession. The higher sedimentation rates are consistent with a greater occurrence of coarse-grained facies and indicated by the thicker stratigraphic package representing a given conodont biozone in the Narao section relative to the Naqing section (Fig. 2). The latter reflects the shallower-water depositional environment of the Narao succession. More precise correlation between the records is not possible given the uncertainty of estimated sedimentation rates, particularly for successions composed of different lithofacies and lacking precise age control.

Although it is possible that the absolute isotopic values of our samples were systematically shifted to more positive values during early diagenesis (i.e., syndepositional carbonate dissolution and recrystallization) (47, 48), we posit that the negative δ238U excursion most likely preserves the original seawater isotopic variability. The Naqing and Narao successions were deposited on a carbonate slope that faced an open-water seaway that was situated at the confluence of the Paleo-Tethys (west) and Panthalassic oceans (east) (). Thus, the depositional environments were exposed to open marine waters with a global marine signal, as supported by the δ13Ccarb record (51). The high δ238U values both below and above the KGB (a range overlapping that of modern shallow-water carbonate sediments) (47, 48) suggest an overall well-oxygenated Carboniferous ocean consistent with high atmospheric pO2 at that time, based on modeling efforts (52–54). Conversely, the δ238U negative excursion across the KGB indicates a shift to more widespread anoxic conditions (Fig. 2).

To constrain the extent of anoxic seafloor area, we modified a community standard U isotope-mass balance model (37, 38) using a stochastic approach for error propagation and explicitly scaled the riverine input of U to increases in the extent of silicate weathering derived from the LOSCAR model (Fig. 2). Specifically, for error propagation, we conducted a Monte Carlo simulation in which parameters in the U isotope-mass balance had uniform distributions reflecting previously proposed values (). Based on these modeling results, we found a substantial increase to ∼18.6% (with 1 SD of 6.6%) in the extent of area of anoxic seafloor is needed to explain the observed negative shift in δ238U (Fig. 3 and ). The model results are consistent with qualitative evidence for a KGB anoxic event inferred from other sedimentary and biotic archives. Widespread occurrence of black shales across the KGB interval in the West Texas Midland Basin (Finis Shale), Midcontinent shelf (Heebner Shale of the Oread Cyclothem), Illinois Basins (middle shale of the Shumway cyclothem), and the Appalachian Basin (Ames Limestone and associated black shale), in the United States (21, 55), as well as the Qian-Gui Basin of South China (in the present study), supports deposition under anoxic conditions. Furthermore, there is a significant (∼25%) drop in biodiversity of benthic faunas (e.g., foraminifers and brachiopods) beginning in the late Moscovian and reaching the nadir immediately below the KGB, which is superimposed on a long-term, late Carboniferous–early Permian biodiversification event (Fig. 1 and ) (56, 57).

For direct comparison, we carried out the same linked LOSCAR–U modeling used for the KGB event for four greenhouse C-perturbation events (end-Permian, end-Triassic, Cretaceous Cenomanian–Turonian anoxic event [OAE2], and PETM) (). There is a relationship between the rate of C release, or the temperature increase, and the extent of anoxia during warming events (Fig. 4). The KGB falls off this trend, with extensive anoxia despite a relatively slow C-release rate (Fig. 4). In other words, we found that the estimated increase in the extent of anoxic conditions during the KGB—relative to the C-release rate or the change in temperature during the perturbation—was more extreme than warming events that occurred under greenhouse conditions. This striking difference (Fig. 4) suggests that the KGB event was primed to develop anoxia and raises the question as to whether there are fundamental differences in redox regulation and in the potential to develop anoxia under a glacial climate state (icehouse) versus for an ice-free greenhouse state.

Fig. 4.

Comparison between the KGB warming event and greenhouse C-perturbation events over the last 300 Myr, indicating that the KGB warming had a more extreme impact on the marine redox landscape relative to its rate of C injection and SST increase than the other events. (A) Increase in fanox versus rate of C injection. (B) Increase in fanox versus rate of temperature (T) rise. The dashed trend lines are plotted based on all events excluding the KGB and PETM. Data points indicate the average of simulations from this study (linked LOSCAR and U modeling) and from previous studies (). Error bars are 1 SD. TJB, Triassic–Jurassic boundary; TOAE, Toarcian ocean anoxic event.

Mechanism of Marine Anoxia.

To interrogate the potential mechanisms that may have been responsible for significant ocean anoxia during the KGB warming, we carried out climate model simulations using the fully coupled Community Earth System Model (CESM), version 1.2 (58), with a KGB event paleogeography. Our simulations suggest that ocean deoxygenation was tied, at least in part, to enhanced thermocline stratification and weakened deep meridional overturning circulation in the Northern Hemisphere. In a glacial state (low CO2 [LowCO2] simulation, 280 ppm), deep-water formation occurs in both hemispheres of the Panthalassic Ocean (Fig. 5 ). Ocean warming in an interglacial state (high CO2 [HighCO2] simulation, 560 ppm) decreases surface seawater density and induces a poleward migration of isopycnal outcrops, as compared with the glacial state, indicating enhanced seawater stratification during warming within the late Paleozoic icehouse (Fig. 5 ). Critically, increased surface stratification reduces mixed-layer depths (Fig. 5 and ), which will decrease ocean-interior oxygen concentrations (Fig. 5) (59). Warming-induced surface stratification in the northern high latitudes also causes a shutdown of deep convection that, in turn, leads to a 61% reduction in maximum overturning in the Northern Panthalassic Ocean (). These results suggest that a restructuring of circulation regimes may have played a role in inducing anoxia, although additional work is necessary to quantify the effects on the global marine redox landscape.

Fig. 5.

Coupled climate model simulations for the KGB. (A and B) The climatological maximum mixed-layer depth, showing decreased northern, high-latitude, deep-water formation and shallower, low-latitude mixed layer in the HighCO2 (560 ppm; interglacial state) compared with LowCO2 (280 ppm; glacial state) simulations. (C and D) Zonal mean density (kg/m3; black-filled contours) and temperature (°C; red contours) in the upper 500 m of the Panthalassic Ocean for the LowCO2 and HighCO2 simulations. Note that isopycnal outcrops are shifted poleward in the HighCO2 simulation as compared with the LowCO2 simulation, suggesting increased seawater stratification with a doubling of pCO2. (E and F) Differences between the LowCO2 and HighCO2 simulations in mean annual sea surface temperature (°C) (E) and in oxygen saturation concentration (mol/m3) (F). Oxygen saturation concentration is calculated according to ref. 80. The empirical formula expresses the dependence of the logarithm of oxygen solubility on temperature and salinity in surface ocean waters. White areas in the high latitudes correspond to SST < 0 °C and thus fall outside the temperature range of the empirical formula.

Additionally, the greater extent of anoxia observed at the KGB per degree of warming, relative to other C-perturbation events that occurred under greenhouse climates, can be attributed to multiple reasons. First, the logarithmic relation between temperature change and atmospheric pCO2 predicts that the same rate of C injection may cause a greater temperature increase when the baseline, prewarming pCO2 is low in icehouse conditions than when prewarming pCO2 is high in greenhouse conditions. Furthermore, the climate sensitivity can also be amplified (by a factor of 2 or more) in icehouse conditions, due to ice-sheet changes that constitute a positive feedback (with stronger impact of sea ice than land ice), relative to greenhouse conditions (60). The greater warming response to pCO2 change under icehouse conditions would lead to more widespread anoxia through reduced O2 solubility and/or weakened ocean ventilation (Fig. 5) (1). Second, abrupt warming during icehouse conditions could have led to a pulse of phosphorus input to the ocean released during weathering of previously glaciated catchments (61, 62), promoting increased primary productivity that would have ultimately elevated oxygen consumption from the water column, enhancing ocean anoxia (via expansion of OMZs). Last, the late Paleozoic continental configuration, and thus the circulation patterns in the Paleo-Tethys Ocean, might have created a large-scale nutrient trap that promoted widespread euxinia in the region (e.g., South China) (63).

Collectively, our findings from the late Paleozoic icehouse provide a unique insight into global warming–induced ocean deoxygenation in a glacial state. More broadly, this work highlights that the extent of deoxygenation, relative to temperature change, can be highly variable for different C-injection events, particularly between those that occurred under different climate states (Fig. 4). Although the drivers of this difference in redox stability in response to warming under different climate states are not well constrained, our study provides a clear motivation for further study of this issue.

Materials and Methods

Carbonate Uranium Isotopes.

Fresh rock samples were cut into small pieces (∼2 × 2 cm) avoiding any altered parts (e.g., calcite veins, chert nodules, weathered surfaces). Samples were pulverized to ∼200 mesh. Minimally 3 g of sample powder was fully digested in 1 M HNO3 in a 50-mL trace-metal-clean centrifuge tube. The resulting supernatant was spiked using a double-spike solution containing 236U and 233U. The spike U was purified using a uranium and tetravalents actinides resin (Eichrom) following procedures described in ref. 64. The δ238U values were measured with a ThermoFinnigan Neptune multicollector–inductively coupled plasma mass spectrometer (ICPMS) instrument at the W. M. Keck Laboratory for Environmental Biogeochemistry, at Arizona State University, and the Yale University Metal Geochemistry Center. The U isotopic compositions of samples are reported relative to those of CRM145 standard. The analytical precision of this method is better than ±0.08‰ (2 SD), based on long-term measurements of the standard CRM145.

Major and Trace Elements.

Major and trace element concentrations were measured on splits of each supernatant using a Thermo iCAP quadrupole ICPMS at the W. M. Keck Laboratory for Environmental Biogeochemistry, Arizona State University. Typical precision was better than 3% and 5% (2 SD) for major and trace elements, respectively, based on repeated analysis of in-run check standards.

Organic C Isotopes.

Approximately 5 g of carbonate powder was fully dissolved in 10% HCl. Undissolved residues were filtered using 45-µm cellulose nitrate filters under vacuum and dried at 40 °C in a laboratory oven. Samples of acid insoluble residues (50 to 120 mg) were loaded in tin boats and analyzed for δ13Corg using an Elementar Vario EL Cube elemental analyzer (Elementar Analysensysteme GmbH, Hanau, Germany) interfaced to an Isoprime VisION IRMS at the University of California, Davis, and a Flash 2000 Elemental Analyzer interfaced to a Delta V Plus IRMS at Washington University, St. Louis, MO. During analysis, samples are interspersed with several replicates of at least four different laboratory reference materials. A sample’s provisional isotope ratio is measured relative to a reference gas peak analyzed with each sample. Provisional values are corrected using the known values of the included laboratory reference materials. The long-term precision for δ13Corg is 0.2‰ (±2 SD). The corrected delta values are expressed relative to the international standard, Vienna Pee Dee Belemnite.

LOSCAR Modeling.

C-cycle simulations for C emissions during the Late Pennsylvanian were performed using the LOSCAR model (31). This model is designed to calculate the partitioning of C among ocean, atmosphere, and sediments on various time scales and is suitable for simulating the atmospheric pCO2 level and the ocean C isotopes during and after C injection over thousands to millions of years. Although the LOSCAR model is not specifically designed to simulate the KGB C-cycle perturbation, it can be modified to better match the specific boundary conditions. We used default modern settings (icehouse), instead of the PETM settings (greenhouse), of the LOSCAR model for the simulation. We made several modifications to the LOSCAR model to (1) match the baseline δ13Ccarb (∼3.8‰) and pCO2 (∼350 ppm) before the KGB event, and (2) accommodate the significantly larger partitioning of shelf CaCO3 rain to shallow-water sediments versus pelagic CaCO3 rain to deep-sea sediments in the late Paleozoic era relative to the modern era.

First, we changed the δ13C value of weathered carbonate on land from 1.5 to 4.5‰, which is similar to the δ13Ccarb used for the time interval straddling the KGB interval in the GEOCARBSULF model (65). Second, we set the baseline pCO2 to 350 ppm and the LOSCAR model was spun up until it reached this baseline value. Third, we changed the ratio of shelf-carbonate precipitate to deep-sea rain from 1 to 5 to reflect the evolution of pelagic calcifiers, with the value within a relevant range (31). We further incorporated in the LOSCAR modeling the conveyor transport value (23.9 Sv prior to the KGB warming and 22.5 Sv during the KGB warming) defined by our CESM, version 1.2, simulation (see Coupled Climate Model Simulations below). To match the timing of initial decrease of measured δ13Ccarb and the initial increase of pCO2, our model ran with a C injection beginning at 304.12 Ma. As the amount and duration of the C emission and the δ13C of the released C are unknown, we combined a grid search with a manual fitting to search for the best estimates of these parameters. The different model configurations also serve as sensitivity tests ().

The emission total was selected from 1,000, 2,500, 5,000, 10,000, 15,000, and 20,000 Gt C. Emission duration was selected from 100, 200, 300, 400, 500, and 600 kyr. The δ13C of the C source was selected from 0, –5, –40, and –60‰, which corresponds to the δ13C value of inorganic thermogenic C, volcanic degassed C, thermogenic organic C, and biogenic C, respectively. In addition, we conducted a grid search for a climate sensitivity of 1.5, 3, and 6. In total, 432 scenarios were used. Model runs were implemented by bridging the open-source language R and the LOSCAR model. By comparing each model run with the observed δ13Ccarb values and CO2 estimates (see for δ13C and pCO2, respectively) using visual inspection and the R2 score of the modeled versus proxy record values for both pCO2 and δ13C, we concluded that the runs with emission totals of 5,000 to10,000 Gt C, emission durations of 200 to 400 kyr, and δ13C values of –5 to –40‰ generally span our pCO2 and δ13C records. Lower C-emission rates (e.g., 1,000 and 2,500 Gt C) produced modeled pCO2 that is too low (), whereas higher C-emission rates (e.g., 15,000 and 20,000 Gt C) resulted in modeled pCO2 that is too high as compared with the pCO2 estimates (). Following the grid search, we further manually varied the C-emission amount, the δ13C of the released C, and the C-emission duration in order to find the best range of these parameters that yields the best scenario through visual inspection and R2 score calculation. This was found to be an emission total of 9,000 Gt C, emission duration of 300 kyr, and δ13C of –19‰ (Fig. 3).

Uranium Cycle Modeling and Estimates of Anoxic Seafloor Area.

We adopted a dynamic U isotope-mass balance model and the measured U-isotope values [mean of the locally weighted scatterplot smoothing (LOWESS) regression subtracting 0.27‰, updated diagenetic correction (48)] from the Naqing section to constrain the extent of seafloor area that was overlain by anoxic bottom water (fanox). Specifically, using a Monte Carlo sampling approach, we solve for fanox and variations to the size of the marine uranium inventory (Nsw) according to the following expression in refs. 37 and 38:which is derived from the following U mass and isotope balance equations:

where

Here, δ238Usw, δ238Uriv, δ238Uanox, and δ238Uother are the U-isotope compositions of seawater, the riverine source (Jriv), anoxic sedimentary sinks (Janox), and remaining other sinks (Jother), respectively. The Jriv of U is scaled linearly with the weathering flux (total silicate and carbonate) based on LOSCAR modeling outputs. Δanox and Δother are the effective fractionation factors associated with anoxic and nonanoxic sediment deposition. The terms kanox and kother are the rate constants and describe U burial rates in average anoxic and all other settings (). Notably, there are differences in kanox values used in U modeling (), which has led to variable magnitudes of fanox in the literature () (38, 66). The differences are mainly derived from the assumption of different values of modern ocean fanox and scavenging rates of U (38, 45, 66, 67). In this study, we ran the model with the more conservative kanox value (5.10E-05 y−1 ± 20%) for all events, and applied the two kanox end-member values (5.10E-05 y−1 and 1.74E-04 y−1) to a sensitivity analysis using a full kanox range (). We ran numerous simulations (n = 10,000), randomly selecting input parameters from a prescribed range (). At the start of each simulation, a spin-up is performed to achieve steady state (i.e., change in Nsw less than 1% over 0.5 Myr) with the selected parameters. Final steady-state values are subsequently selected as inputs for simulation of the event.

For greenhouse OAEs with available carbonate δ238U data [e.g., the Permian–Triassic boundary (37, 41, 64), Triassic–Jurassic boundary (68), OAE2 (16), as well as the PETM (49)], we ran the same linked LOSCAR and U models in two scenarios of low and high C injection as reported in the literature (). Prior to running the LOSCAR model, the boundary conditions for these greenhouse events were modified to match their initial pCO2 levels prior to the C-isotope excursion. For the Jurassic (Toarcian) OAE, we utilized information about the anoxic seafloor area (∼5%) based on molybdenum isotopes (69), adding a 50% error for uncertainty. Finally, we plotted changes in fanox against the rate of C injection and rate of temperature rise for all of the events (Fig. 4 and ) in order to evaluate the impact of C injection and warming on the marine redox landscape under background icehouse versus greenhouse climate states.

Coupled Climate Model Simulations.

Model description.

The simulations were performed with the fully coupled CESM, version 1.2. The fully coupled configuration of CESM includes the Community Atmosphere Model v(CAM), version 5; Community Land Model (CLM), version 4; Parallel Ocean Model (POP), version 2; and Sea Ice Model (CICE), version 4, components (58). CAM and CLM have a 1.9° × 2.5° horizontal resolution with 30 vertical levels. POP and CICE have a nominal 1° horizontal resolution. POP has 60 vertical levels with 10-m resolution in the top 15 levels that increase to 250-m resolution for the deepest ocean. Vertical mixing was represented using the modified K-Profile Parameterization scheme (70) with a horizontally constant background diffusivity of 0.1 × 10−4 m2 s−1 to accommodate the deep-time paleogeography. The Gent–McWilliams scheme with diagnostically evaluated isopycnal diffusivity, as described in ref. 71, was used as the lateral closure with upper-ocean isopycnal diffusivity values up to 3,000 m2 s−1. The horizontal diffusivity coefficient was also set to 3,000 m2 s−1 in the surface diabatic layer. A detailed description of the ocean model can be found in ref. 71.

Experimental design.

We performed two Earth system model simulations of late Pennsylvanian, referred to as LowCO2 and HighCO2, which were each run for ∼2,500 y. At the end of the integration, the simulations reached quasi-equilibrium in the deep ocean, as characterized by time series of temperature and meridional overturning circulation (). The LowCO2 and HighCO2 simulations include the same continental distribution of a late Pennsylvanian (300 Ma) paleogeography (72), deep-ocean bathymetry, treatment of aerosols, and solar luminosity. Deep-ocean bathymetry includes an idealized midocean ridge system defined by accretion zones in the Gplates tectonic reconstruction at 300 Ma (http://www.gplates.org). Midocean ridge crests have a uniform depth of 2,500 m, based on the present-day global average (73). Elsewhere the abyssal seafloor has a uniform depth of 4,000 m and continental shelves at a depth of 60-m. Aerosol distributions were zonally averaged separately for land and ocean from preindustrial levels following the procedure outlined by Heavens et al. (74). We used a modern Earth orbital configuration, and the total incoming solar irradiance was reduced to a late Pennsylvanian value of 1,333 W m−2 (97.5% of present), following Gough (75).

The late Pennsylvanian LowCO2 and HighCO2 simulations differ with respect to atmospheric pCO2, tropical vegetation, ice sheets, and sea level to represent key differences between interglacial highstand and glacial lowstand conditions during the late Pennsylvanian. LowCO2 and HighCO2 atmospheric CO2 concentrations were specified as 280 ppm and 560 ppm, respectively, based on a multiproxy (103- to 104-y resolution) CO2 reconstruction (21, 24). In the absence of proxy data, other greenhouse gas concentrations (e.g., N2O, CH4) were set to their preindustrial values. The distribution of prescribed vegetation follows the mid to late Pennsylvanian biomes for interglacial and glacial phases (76), where the closest possible composition of modern-plant functional types was chosen to represent each biome. The prescribed land ice geographies for LowCO2 and HighCO2 simulations are based on reconstructions of paleo-glacier extent and ice-stream pathways constrained by Carboniferous glacial features and deposits from high-latitude southern Gondwana as well as their detrital zircon U–Pb geochronology (77, 78). In our simulations, continental ice was prescribed by lowering the albedo of the land surface, and the extent of elevation of the ice sheets is fixed in time. The CESM, version 1.2, we used is not biogeochemistry enabled; thus, dissolved oxygen was not included in the LowCO2 and HighCO2 simulations. In this study, we chose to model the mean annual oxygen-saturation concentration difference between the HighCO2 and LowCO2 simulations (Fig. 5) rather than explicitly simulating dissolved oxygen, because there are significant uncertainties that arise when applying the modern plankton functional types (used in the biogeochemistry-enabled CESM) and the associated parameterizations of primary productivity to the KGB event, which predated the evolution of modern plankton.