• Cretaceous/Paleogene boundary
• ocean acidification
• boron isotopes
• mass extinction
• GENIE model

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 (SI Appendix). Importantly, one shallow marine sample site studied here (Geulhemmerberg Cave) is thought to preserve sediments from the first 10² to 10³ 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 δ¹¹B (δ¹¹B_sw; see SI Appendix), one of the main uncertainties in calculating pH from deep-time δ¹¹B measurements. Assuming feasible limits on the calcium carbonate saturation state (Ω_CaCO3) and oxygen utilization at each site (see SI Appendix, section 5C), we constrain K-Pg δ¹¹B_sw 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. 1C). These data provide geochemical evidence of rapid acidification in the immediate aftermath of the Chicxulub impact. Combined with model-derived estimates of alkalinity (SI Appendix), this change in pH corresponds to a rise in atmospheric partial pressure of CO₂ (pCO₂) from ∼900 ppm in the latest Maastrichtian to ∼1,600 ppm in the immediate aftermath of bolide impact (Fig. 1D). 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 SO₂ and NO_x 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 CO₂ release would have been required. This additional CO₂ 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 SI Appendix), and/or a contribution from a pulse of volcanism (5).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

Records of surface ocean foraminiferal δ¹¹B (B), calculated pH (C), and pCO₂ (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 δ¹¹B_sw, 39.45 ± 0.4‰. pCO₂ 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 pCO₂ can be seen in Fig. 2 and in Dataset S1. For details of the age models used, the estimation of δ¹¹B_sw, carbonate system calculations, and uncertainty propagation see SI Appendix.

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 δ¹³C values of surface- and deep-dwelling calcifiers (e.g., ref. 23). Paired measurements of δ¹¹B (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 δ¹³C 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 δ¹³C 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 δ¹³C 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 δ¹³C 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 pCO₂ (thousands of parts per million) could have driven convergence in deep- and surface-ocean δ¹³C 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. 1A). 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 CO₂ degassing in response to impact (4), its effects were overwhelmed by the biogeochemical effect of extinction. Other independent pCO₂ proxy records (e.g., ref. 33) likewise do not support a substantial rise in CO₂, and there is no evidence for deep-sea warming at this time (34). Thus, we discount high pCO₂ as the driver of reduced vertical δ¹³C gradients (29) and infer that the convergence in surface and deep-ocean δ¹³C 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 δ¹³C gradients, we used GENIE (Grid ENabled Integrated Earth System Model) (35), with a late Maastrichtian bathymetry and continental configuration (SI Appendix, Fig. S17 and SI Appendix). We simulate the response of marine δ¹³C 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 SI Appendix). For each simulation, we show depth transects for δ¹³C and pH at Shatsky Rise (ODP Site 1209) (Fig. 2 B and C). Pacific basin δ¹³C 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 δ¹³C gradient with surface ocean δ¹³C values >1‰ lighter than benthic (Fig. 2B) and a muted but positive surface-deep pH gradient (Fig. 2C). The “Living Ocean” scenario with shallower organic matter remineralization leaves a muted but positive surface-deep δ¹³C gradient, but with a flattened pH gradient. Only with a 50% reduction in global new primary productivity can our model produce a flattened vertical δ¹³C gradient together with a positive pH gradient.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Paired benthic–planktic foraminiferal δ¹¹B 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 δ¹¹B/pH in the surface ocean relative to the deep until 65.92 Ma, where this gradient disappears (see also SI Appendix, Fig. S12). Shown for context is the convergence in δ¹³C from surface waters (bulk CaCO₃, in green) and deep ocean waters (benthic foraminifera, navy) (A). Water column δ¹³C 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 δ¹³C 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 SI Appendix for more details.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Pacific zonal means of the simulated depth distribution of δ¹³C_DIC 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 CO₂ at high latitudes due to increased solubility of surface waters with respect to CO₂ as they cool, imparts an inverse vertical δ¹³C_DIC 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 δ¹³C gradient (Fig. 2B) rather than the observed approximately flat gradient. This phenomenon arises because of the solubility pump, that is, the increased solubility of CO₂ (and hence CO₂ uptake) in cooling surface waters as they reach high latitudes, where they ultimately sink to form deep water. Because the fractionation between dissolved CO_2(aq) (and thus atmospheric pCO₂) and aqueous bicarbonate (and thus dissolved inorganic carbon, DIC) increases by about 0.1‰ per °C cooling, the surface δ¹³C_DIC at colder, higher latitudes is driven higher with respect to the atmosphere compared to the warmer, low-latitude ocean surface (Fig. 3C). 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 δ¹³C, with most positive values at the surface (Fig. 3A and SI Appendix). 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 δ¹³C gradient opposite in sign and equal in magnitude to the solubility pump, canceling out the difference in δ¹³C 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 δ¹³C) lens of southern-sourced water (Figs. 2B and 3B). 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 δ¹³C gradients at Shatsky Rise (Fig. 2A) 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 δ¹³C gradients without an associated collapse in pH gradients (Fig. 2A) are consistent with a substantial reduction in new primary productivity globally following mass extinction (Fig. 2 B and C and SI Appendix). The magnitude of productivity loss required depends on changes in carbon isotope vital effects in pelagic CaCO₃ producers, which likely account for some (but not all) of the observed collapse in δ¹³C 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 (SI Appendix, Figs. S23 and S24).

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. 2A). 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. 2C). Shallow water remineralization alone (i.e., a “Living Ocean” scenario) does not fully collapse the δ¹³C gradient in our model, but combined with known carbon-isotope vital effects that drive surface-derived carbonates lighter than ambient δ¹³C_DIC (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 δ¹³C and pH gradients between the surface and deep sea.

Reduced vertical δ¹³C_CaCO3 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 δ¹¹B data do not support prolonged surface ocean acidification, making it difficult to attribute this delayed biogeochemical recovery to acidification from sustained (or enhanced) CO₂ 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.

• Copyright © 2019 the Author(s). Published by PNAS.