Rapid ocean acidification and protracted Earth system recovery followed the end-Cretaceous Chicxulub impact.

Mass extinction at the Cretaceous-Paleogene (K-Pg) boundary coincides with the Chicxulub bolide impact and also falls within the broader time frame of Deccan trap emplacement. Critically, though, empirical evidence as to how either of these factors could have driven observed extinction patterns and carbon cycle perturbations is still lacking. Here, using boron isotopes in foraminifera, we document a geologically rapid surface-ocean pH drop following the Chicxulub impact, supporting impact-induced ocean acidification as a mechanism for ecological collapse in the marine realm. Subsequently, surface water pH rebounded sharply with the extinction of marine calcifiers and the associated imbalance in the global carbon cycle. Our reconstructed water-column pH gradients, combined with Earth system modeling, indicate that a partial ∼50% reduction in global marine primary productivity is sufficient to explain observed marine carbon isotope patterns at the K-Pg, due to the underlying action of the solubility pump. While primary productivity recovered within a few tens of thousands of years, inefficiency in carbon export to the deep sea lasted much longer. This phased recovery scenario reconciles competing hypotheses previously put forward to explain the K-Pg carbon isotope records, and explains both spatially variable patterns of change in marine productivity across the event and a lack of extinction at the deep sea floor. In sum, we provide insights into the drivers of the last mass extinction, the recovery of marine carbon cycling in a postextinction world, and the way in which marine life imprints its isotopic signal onto the geological record.

There is abundant evidence for a massive bolide impact at the Cretaceous–Paleogene (K-Pg) boundary (66.04 Ma) that coincides in those sections with sufficient temporal resolution with a mass extinction horizon (e.g., ref. 1). Geological records document many effects from impact, including massive tsunamis, earthquake-driven gravity flows, and molten ejecta fallout (e.g., ref. 1). However, the mechanism (or mechanisms) by which impact drove global-scale ecosystem turnover and mass extinction is less certain. Among the most prominent hypotheses are global darkness and associated primary productivity loss leading to food chain collapse, acid rain, impact winter, and flash ocean acidification (2). Some of these mechanisms are supported by modeling work (e.g., ref. 3) but, critically, they generally lack empirical validation. Furthermore, the issue is complicated by the possibility of contributing effects from ongoing or intensified Deccan flood basalt volcanism (4, 5) and the short timescales of at least some of the effects of impact compared to the resolution of the geological record.

Here we apply the boron isotope-pH proxy to planktic and benthic foraminifera to examine whether ocean acidification occurred at the K-Pg boundary and to test competing hypotheses that have been proposed to explain changes in the marine carbon cycle in the aftermath of the mass extinction. We analyzed material from multiple sites where sufficient well-preserved foraminifera from the mixed layer of the surface ocean could be obtained, with >7,000 small, thin-walled postextinction foraminifera required per analysis in some cases. By combining measurements from continental shelf sites in communication with the open ocean [Geulhemmerberg Cave, The Netherlands (6), Owl Creek, Mississippi (7), and Brazos River, Texas (8)] and deep-sea drilling sites (Deep-Sea Drilling Program Site 465 and Ocean Drilling Program [ODP] Site 1209 in the Pacific and ODP Site 1049 in the Atlantic), we construct a composite global record spanning ∼800 ky (). Importantly, one shallow marine sample site studied here (Geulhemmerberg Cave) is thought to preserve sediments from the first 102 to 103 y after bolide impact (6), allowing us to pinpoint a short-lived signal of impact-induced ocean acidification that cannot be temporally resolved in deep marine records. Samples from ODP Site 1209 allow us to investigate longer-term changes in the ocean’s biological carbon pump as expressed in vertical pH gradients, as this site is one of few deep marine K-Pg boundary sites that preserves sufficient numbers of well-preserved epifaunal benthic and surface-dwelling planktic foraminifera for boron isotope analysis. The broad range of in situ temperatures, hydrostatic pressures, and environmental conditions represented in our sample set allows us to constrain bulk seawater δ11B (δ11Bsw; see ), one of the main uncertainties in calculating pH from deep-time δ11B measurements. Assuming feasible limits on the calcium carbonate saturation state (ΩCaCO3) and oxygen utilization at each site (see ), we constrain K-Pg δ11Bsw to within 39.05 to 39.85‰.

Our data suggest ocean pH was stable over the last 100 ky of the Cretaceous, despite suggestions of a major pulse of Deccan volcanism within this time interval (5). However, our data from the Geulhemmerberg Cave boundary clay indicate a marked ∼0.25 pH unit surface ocean acidification event within a thousand years (6) of the bolide impact in the Geulhemmerberg Cave boundary clay (Fig. 1). These data provide geochemical evidence of rapid acidification in the immediate aftermath of the Chicxulub impact. Combined with model-derived estimates of alkalinity (), this change in pH corresponds to a rise in atmospheric partial pressure of CO2 (pCO2) from ∼900 ppm in the latest Maastrichtian to ∼1,600 ppm in the immediate aftermath of bolide impact (Fig. 1). This measured pH drop may be a conservative estimate of the full magnitude of impact-induced acidification, as the earliest Ir-rich ejecta layer at Geulhemmerberg is thought to have been eroded within ∼100 y of the impact (6). A geologically rapid 0.2 to 0.3 pH unit change would have disadvantaged calcifying plankton vs. noncalcifiers (9) and could therefore explain the selective extinction of calcifying pelagic organisms during the K-Pg mass extinction (10–13) compared to noncalcareous groups such as dinoflagellates (e.g., ref. 14) and radiolarians (e.g., ref. 15). Similarly, since coccolithophores and planktic foraminifera from coastal areas may be better adapted to fluctuations in pH (16), and smaller foraminifera may better maintain calcification under low pH (17), acidification could also help to explain on-shelf/off-shelf and size-selective patterns of extinction among calcifiers (18, 19). The observed magnitude of acidification is compatible with the rainout of SO2 and NOx released by the impactor (10, 20, 21), depending on the timing of sediment deposition at Geulhemmerberg. The effect of rainout products was transient (e.g., ref. 21), and if Geulhemmerberg clays were deposited on timescales greater than a few hundred years, a contribution to acidification from CO2 release would have been required. This additional CO2 could be sourced from carbonate target rocks (ref. 1 and references therein), the biosphere [through wildfires and decay (ref. 1 and references therein)], a transient weakening of the biological carbon pump (ref. 22 and references therein and ), and/or a contribution from a pulse of volcanism (5).

Fig. 1.

Records of surface ocean foraminiferal δ11B (B), calculated pH (C), and pCO2 (D) across the K-Pg boundary, with high-resolution foraminiferal diversity counts from ref. 39 at the K-Pg Global Boundary Stratotype Section and Point (El Kef) plotted for context (A). pH is calculated assuming our best estimate of K-Pg δ11Bsw, 39.45 ± 0.4‰. pCO2 is calculated from pH along with total alkalinity estimates at each site from a GENIE late Maastrichtian simulation, adjusted for dynamic changes in alkalinity across the K-Pg using LOSCAR simulations from ref. 22 that match observed patterns of carbonate burial. Gray shaded areas are 1-sigma uncertainties, with thin lines representing 1,000 Monte Carlo simulations from the 10,000 that were run. For clarity, we only plot those samples that represent the surface mixed layer, which should be approximately in equilibrium with the atmosphere. Additional data from deep-sea benthic and thermocline-dwelling planktic foraminifera that do not reflect atmospheric pCO2 can be seen in Fig. 2 and in Dataset S1. For details of the age models used, the estimation of δ11Bsw, carbonate system calculations, and uncertainty propagation see .

Fig. 2.

Paired benthic–planktic foraminiferal δ11B measurements from 4 time slices following the K-Pg boundary at ODP Site 1209 (Shatsky Rise) are used to calculate the difference in surface and deep ocean pH (∆pH; red diamonds, A). These data show a trend from “normal” patterns of higher δ11B/pH in the surface ocean relative to the deep until 65.92 Ma, where this gradient disappears (see also ). Shown for context is the convergence in δ13C from surface waters (bulk CaCO3, in green) and deep ocean waters (benthic foraminifera, navy) (A). Water column δ13C and pH profiles in the cGENIE model the location of Shatsky Rise (see also Fig. 3) are shown in B and C, respectively. A weakened or reversed δ13C gradient can arise in a number of ways (B), but only a shallow remineralization (or “Living Ocean,” red dot-dashed line) model can produce a flat gradient in pH between the surface 80 m and 2,500-m water depth (C). Carbonate data from ODP 1209 include new data and benthic data from ref. 34. See for more details.

Loss of plankton functional groups with the K-Pg impact has profound implications for marine carbon cycling (22), manifest in one of the earliest and most striking geochemical observations of the K-Pg boundary in open ocean sediments: a convergence in the δ13C values of surface- and deep-dwelling calcifiers (e.g., ref. 23). Paired measurements of δ11B (and hence pH) from surface and deep water after the K-Pg boundary allow us to address a long-running, unresolved debate about the cause of this convergence, and the long-term effect of the K-Pg mass extinction on marine carbon cycling. Surface ocean δ13C is typically maintained at heavy values relative to the deep sea through the biological pump, which exports isotopically light organic carbon from the surface ocean to depth. The sharp drop in surface ocean δ13C at the K-Pg was therefore initially interpreted as a near-total loss of marine new primary productivity after impact, in the “Strangelove Ocean” hypothesis (24). Later, the “Living Ocean” hypothesis (25) posited that such a reduction of vertical δ13C gradients arose instead from a decrease in the efficiency of the biological pump, that is, continued productivity but reduced sinking efficiency leading to more remineralizaton of organic carbon in the upper water column. Evidence for export to the deep ocean persisting in at least some regions led to the “Heterogeneous Ocean” hypothesis (11, 26), while it has also been argued that part of the reduction in δ13C gradients stems from the fact that carbon isotope vital effects differ in postextinction carbonate producers compared to the preextinction ones (11, 27). These latter scenarios appear more compatible with low levels of benthic foraminiferal extinction at the K-Pg (11, 28) than a near-total loss of primary or export production. Another alternative hypothesis not linked to any marine biogeochemical response is that extremely high atmospheric pCO2 (thousands of parts per million) could have driven convergence in deep- and surface-ocean δ13C inorganically (29). Such a scenario might arise if outgassing from the Deccan Traps was accelerated following the impact (4).

Our records show that initial surface acidification was followed by a rapid rebound and overshoot in surface ocean pH (∆pH ≈ 0.5) within 40 ky of the K-Pg boundary, returning to preevent values after an additional ∼80 ky (Fig. 1). This is consistent with model predictions of ocean pH rise due to the extinction of calcifying organisms and a consequent transient reduction in the marine alkalinity sink (22, 30). Such a period of imbalance in marine alkalinity fluxes in the earliest Danian is supported by observations of improved deep-sea carbonate preservation (22), the occurrence of heavily calcified deep-sea benthic foraminifera (31), and the deposition of abiogenic calcite crystals in the oceans (32). The absence of any prolonged (>50 ky) surface ocean acidification after the K-Pg suggests, therefore, that if there were a state shift toward significantly more active Deccan CO2 degassing in response to impact (4), its effects were overwhelmed by the biogeochemical effect of extinction. Other independent pCO2 proxy records (e.g., ref. 33) likewise do not support a substantial rise in CO2, and there is no evidence for deep-sea warming at this time (34). Thus, we discount high pCO2 as the driver of reduced vertical δ13C gradients (29) and infer that the convergence in surface and deep-ocean δ13C arose from a change in oceanic biogeochemical cycling, combined to some extent with changes in carbonate vital effects.

To explore how changes in the carbon cycle would manifest in vertical pH and δ13C gradients, we used GENIE (Grid ENabled Integrated Earth System Model) (35), with a late Maastrichtian bathymetry and continental configuration ( and ). We simulate the response of marine δ13C and pH given either reduced or zero primary productivity (a “Strangelove Ocean”-type scenario) or less efficient biological pump (i.e., enhanced remineralization in the upper water column, “Living Ocean” scenario; full details are given in ). For each simulation, we show depth transects for δ13C and pH at Shatsky Rise (ODP Site 1209) (Fig. 2 ). Pacific basin δ13C zonal means for 3 illustrative states of biological pumping are shown in Fig. 3. We use the observed patterns of boron and carbon isotope change to distinguish between the competing hypotheses for postextinction biogeochemical change. The fingerprint of the “Strangelove” scenario with a complete loss of productivity is a reversed vertical δ13C gradient with surface ocean δ13C values >1‰ lighter than benthic (Fig. 2) and a muted but positive surface-deep pH gradient (Fig. 2). The “Living Ocean” scenario with shallower organic matter remineralization leaves a muted but positive surface-deep δ13C gradient, but with a flattened pH gradient. Only with a 50% reduction in global new primary productivity can our model produce a flattened vertical δ13C gradient together with a positive pH gradient.

Fig. 3.

Pacific zonal means of the simulated depth distribution of δ13CDIC in the Maastrichtian ocean for (A) the default (full-strength) biological pump, (B) global export production reduced to 50% of the default, and (C) no biological pump. The approximate paleo-latitude of ODP Site 1209 (Shatsky Rise) is marked by the vertical dashed white line. The action of the solubility pump, which absorbs atmospheric CO2 at high latitudes due to increased solubility of surface waters with respect to CO2 as they cool, imparts an inverse vertical δ13CDIC gradient at low latitudes in the absence of biological productivity (C; see the main text for explanation).

The “Strangelove” scenario produces a pronounced >1‰ reversed vertical δ13C gradient (Fig. 2) rather than the observed approximately flat gradient. This phenomenon arises because of the solubility pump, that is, the increased solubility of CO2 (and hence CO2 uptake) in cooling surface waters as they reach high latitudes, where they ultimately sink to form deep water. Because the fractionation between dissolved CO2(aq) (and thus atmospheric pCO2) and aqueous bicarbonate (and thus dissolved inorganic carbon, DIC) increases by about 0.1‰ per °C cooling, the surface δ13CDIC at colder, higher latitudes is driven higher with respect to the atmosphere compared to the warmer, low-latitude ocean surface (Fig. 3). Although this effect of the solubility pump is always there, typically the action of the biological pump overprints the pattern of the solubility pump to give the observed “normal” sign of δ13C, with most positive values at the surface (Fig. 3 and ). Somewhere between a “Strangelove”-style shutdown and a full-strength biological pump, there must exist a weaker state of the biological pump that would impart a δ13C gradient opposite in sign and equal in magnitude to the solubility pump, canceling out the difference in δ13C between surface and deep ocean. We find this occurs with a 50% reduction in global export production for all but the highest-latitude Northern Pacific, where downwelling waters overlay an older (and lower δ13C) lens of southern-sourced water (Figs. 2 and 3). This partially maintained productivity is consistent with the survival of benthic foraminifera at the K-Pg (11, 28).

The combination of our cGENIE simulations and the observed sequence of changes in vertical pH and δ13C gradients at Shatsky Rise (Fig. 2) allow us to propose a hypothesis for the behavior of the marine carbon cycle after the K-Pg impact. During the first tens of thousands of years following impact, reduced vertical δ13C gradients without an associated collapse in pH gradients (Fig. 2) are consistent with a substantial reduction in new primary productivity globally following mass extinction (Fig. 2 and ). The magnitude of productivity loss required depends on changes in carbon isotope vital effects in pelagic CaCO3 producers, which likely account for some (but not all) of the observed collapse in δ13C gradients (11, 27, 36). Reduced global export productivity can, however, align with published observations of a “Heterogenous Ocean” with stable or increased productivity at certain sites [e.g., Shatsky Rise (11, 26)]. A global reduction in new primary productivity may have led to a buildup of nutrients in the surface ocean, and thus some locations (e.g., highly oligotrophic gyres) may have supported increased new/export primary productivity in the Danian relative to their preimpact state, even while export productivity globally was reduced ().

After this period of globally reduced, but not collapsed, export productivity, we see a convergence in shallow and deep-ocean pH in our final time-slice ∼120 ky post-K-Pg (Fig. 2). In our model scenarios, only by remineralizing organic carbon shallower in the water column can one produce similar pH values in both the deep sea and the subsurface (Fig. 2). Shallow water remineralization alone (i.e., a “Living Ocean” scenario) does not fully collapse the δ13C gradient in our model, but combined with known carbon-isotope vital effects that drive surface-derived carbonates lighter than ambient δ13CDIC (27, 36, 37) the expression of such a scenario in the fossil record could be a full collapse. We therefore suggest that after an initial period (up to 40 ky) of globally partially weakened new primary productivity globally, it progressively began to recover, but without full restoration of pelagic ecosystem structure and function. In such a scenario, persistently ineffective ballasting of marine organic matter or other factors could have then resulted in more organic material being remineralized shallower in the water column [i.e., a “Living Ocean” scenario (25)], thereby reducing δ13C and pH gradients between the surface and deep sea.

Reduced vertical δ13CCaCO3 gradients persisted for over 1 My after the K-Pg impact (e.g., ref. 27), which suggests that the efficiency of organic carbon export to the deep sea, a key pelagic ecosystem function, was reestablished slowly. Our δ11B data do not support prolonged surface ocean acidification, making it difficult to attribute this delayed biogeochemical recovery to acidification from sustained (or enhanced) CO2 degassing from the Deccan Traps. Rather, there may be intrinsic constraints on the time required to recover normal marine ecosystem function after such severe global perturbations, despite the short generation times that should make marine plankton ideally suited to rapid evolutionary radiation (38). In this way, the K-Pg event shows that even geologically rapid ocean acidification events can have prolonged and profound biotic ramifications.