Amplifying factors leading to the collapse of primary producers during the Chicxulub impact and Deccan Traps eruptions

The latest Cretaceous (K; Maastrichtian) through earliest Paleogene (Pg; Danian) interval is a time marked by one of five major mass extinctions in the Earth’s history. The synthesis of published data permits the temporal correlation of the K-Pg boundary crisis with two major geological events: 1) the Chicxulub impact discovered in the Yucatán peninsula (Mexico) and 2) the Deccan Traps located on the west-central India plateau. In this paper, environmental and biological consequences from the Chicxulub impact and the Deccan continental flood basalts have been explored using a climate-carbon-biodiversity coupled model called the ECO-GEOCLIM model. The novelty of this study is to investigate how abiotic factors (temperature, pH, and calcite saturation state) act on various marine organisms to determine the primary productivity and biodiversity changes in response to a drastic environmental change. This study shows that the combination of the Deccan volcanism with a 10-km impact led to global warming (3.5°C) caused by rising carbon dioxide (CO2) concentration (+470 ppmv), interrupted by a succession of short-term cooling events, provided by a ‘shielding effect’ due to the formation of sulphate aerosols. The consequences related to these climate changes are the decrease of the surface ocean pH by 0.2 (from 8.0 to 7.8), while the deep ocean pH dropped by 0.4 (from 7.8 to 7.4). Without requiring any additional perturbations, these environmental disturbances led to a drastic decrease of the biomass of calcifying species and their biodiversity by about 80%, while the biodiversity of noncalcifying species was reduced by about 60%. We also suggest that the short-lived acidification caused by the Chicxulub impact, when combined with the Deccan traps, may explain the severity of the extinction among pelagic calcifying species.


Introduction
The Cretaceous-Paleogene (K-Pg) boundary is the best-documented mass extinction event, and calcified pelagic plankton provides the most striking example of biotic loss. Among the 131 Late Maastrichtian species of calcareous nannoplankton (i.e., photosynthetic algae producing minute 1-20 μm calcite platelets), only 10 to 12 survived into the Danian period (Bown, 2005;Bown et al., 2004), corresponding to a turnover of ~93% (Bown, 2005). The planktic foraminifera also underwent an important turnover of 70% to 80% (Macleod et al., 1997;Molina et al., 1998). The severity of extinction among pelagic calcifying species compared to non-calcareous taxa (~50% siliceous diatoms; Macleod et al., 1997) or noncalcifying phytoplankton suggests that a massive acidification of the oceans occurred at the K-Pg boundary (Keller, 2014;Tyrrell et al., 2015). Prinn and Fegley (1987) were the first to advocate the assumption of a wide dispersal of acid rain (HNO 3 ) having a large effect, whereas, using estimates of SO 2 released by the K-Pg impact, D'Hondt et al. (1994) were the first to propose that sulphur dioxide (SO 2 )-induced seawater acidification was an efficient killing mechanism (Sigurdsson et al., 1991(Sigurdsson et al., , 1992Pope et al., 1993;Pierazzo et al., 2003). Ohno et al. (2014) simulated an intense and short-lived acidification to reduce the calcite saturation state of the surface ocean, assuming a rapid scavenging of sulphuric acid aerosols (acid rain). A few months later, Tyrrell et al. (2015) contested this finding, demonstrating that a 10-km impactor was too small to repeat the expected decrease in pH.
The bolide scenario was widely accepted in the 1990s when research advances in radiometric dating of large igneous provinces revealed a close temporal relationship between the emplacement of continental flood basalts and oceanic plateaus with the mass extinctions over the last 380 Ma (Courtillot and Renne, 2003). This close temporal match suggests the existence of a cause-and-effect link between the emplacement of the Deccan Traps and the K -Pg mass extinction (Chenet et al., 2009;Courtillot et al., 1986;Courtillot and Renne, 2003;Wignall, 2001), a result challenging the bolide scenario. With an estimated ~2×10 6 km 3 of extruded lava, it has been 3 35 40 45 50 proposed that Deccan eruptions have injected massive quantities of gases into the atmosphere over a period not exceeding 700 kyrs (Schoene et al., 2015), spanning the K-Pg boundary.
To assess the consequences of the Deccan Traps on the inorganic carbon cycle, Dessert et al. (2001) used a biogeochemical model and were the first to highlight a significant pH decrease in seawater (ΔpH = 0.3) associated with the rise of the partial pressure of carbon dioxide ( pCO 2 ; ΔpCO 2 = +1000 ppmv). The assumptions proposed by Dessert et al. (2001) about the Deccan Traps: 1) the volcanism dynamics (linear and continuous degassing), 2) duration of emplacement (100 or 1000 kyrs), and 3) mass of released gases (70,000 GtCO 2 ) are become obsolete with research updates. In addition, recent findings have demonstrated the importance of SO 2 in the Earth's climate during the Deccan Trap emplacements (Mussard et al., 2014;Schmidt et al., 2016). Considering for the first time the combination of 4090 GtCO 2 and 3200 GtSO 2 , Henehan et al. (2016) found a decline of surface seawater pH close 0.2. They demonstrated that the efficiency of weathering accompanying the induced warming rapidly balances (50 kyrs) the acidification of seawater. Based on this finding, they concluded that the initial environmental perturbation caused by direct gas releases was too weak to explain the major disorders observed in the carbonate record at the K-Pg boundary.
Following the conclusions made by Henehan et al. (2016), this study explores the feedback on inorganic climatic interplays between SO 2 and CO 2 and organic tolerance of marine primary producers to amplify this initial perturbation. Consequently, a climate-carbon model has been developed, including carbon and sulphur cycles coupled with a biodiversity model describing the behaviour of marine ecosystems. Instead of a standard approach where the primary producers are only represented as a single pool of organic carbon driven by the amount of nutrients available in the ocean (Dessert et al., 2001;Henehan et al., 2016;Petersen et al., 2016;Schmidt et al., 2016;Tobin et al., 2017;Tyrrell et al., 2015), here the biodiversity model splits the oceanic biosphere into n species (n = 50, see SI), each characterised by biological factors (such as death and reproduction 4 60 65 70 75 rates) combined with the sensitivity to abiotic factors (seawater temperature, pH, and calcite saturation state). This paper presents the influence of the Chicxulub asteroid and Deccan Traps at the K-Pg boundary using a carbon-climate-biodiversity model (ECO-GEOCLIM). The key points addressed in this study include the quantification of environmental changes provided by the Chicxulub asteroid and Deccan Traps and their respective effects on the marine biodiversity, highlighting the amplifying response to an initial perturbation caused by inorganic and organic carbon cycles.
The paper is structured as follows. Section 2 describes the modelling approach and the novelty compared to previous modelling studies. Section 3 presents the results of the scenario combining Deccan Traps and/or the bolide impact to identify explanations for the marine fauna turnover at the K-Pg boundary. Subsequently, a comparison is made between the experiment findings and available data and the previous modelling results in Section 4. Finally, the model limitations are discussed in Section 5, before a summary of the key findings is presented in Section 6.

Model Description
In this paper, by coupling GEOCLIM and a biodiversity model, a new climate-carbon-biodiversity coupled model (ECO-GEOCLIM) is presented and tested. GEOCLIM (Donnadieu et al., 2006) is a comprehensive climate-carbon model coupling a general circulation model called the fast ocean atmospheric model 1.5 (FOAM 1.5; Jacob, 1997) with a model of the main biogeochemical cycles called COMBINE (Goddéris and Joachimski, 2004). The atmospheric component of FOAM is a parallel version of the National Center for Atmospheric Research's Community Climate Model 2 (CCM2) with upgraded radiative and hydrologic physics incorporated in CCM3 v.3.2 (Jacob, 1997). All simulations have been performed with R15 spectral resolution (4.5° × 7.5°) and 18 vertical levels. In addition, FOAM is used in mixed-layer mode, meaning the atmospheric model is linked 5 85 90 95 100 105 to a 50 m mixed-layer ocean with heat transport parameterised through diffusion. The sea ice module uses the thermodynamic component of the CSM1.4 sea ice model, which is based on the Semtner three-layer thermodynamic snow/ice model (Semtner, 1976). No land ice component was included. It is noteworthy that FOAM does not include an interactive chemistry module; all sulphate reactions and climate forcing are computed in COMBINE (for extended details, see Mussard et al., 2014).
COMBINE is a box model that describes the biogeochemical cycles for carbon, phosphorus, alkalinity, oxygen, and sulphates ( fig.1 electronic supplementary material). The present version was updated to include the seafloor weathering process as an additional long-term carbon sink (Brady and Gíslason, 1997;Coogan and Dosso, 2015;Krissansen-Totton and Catling, 2017). We divided the oceanic crust into six layers from the seafloor to a depth of 500 m in the crust, according to a formalism developed by Charnay et al. (2017). Deep water percolating into the oceanic crust follows the present-day temperature gradient (around 116 K/km), starting from the deep ocean temperature at the water-crust interface. To balance the additional carbon consumption from seafloor weathering, the solid Earth degassing rate was increased by 10% compared to its presentday value (i.e. 6.8×10 12 moles of CO 2 /year for aerial volcanoes and 1.8×10 12 moles of CO 2 /year for mid-ocean ridges) in agreement with the fluxes presented by Wallmann (2001).
For the continental carbon sink, field work shows that silicate weathering is linearly dependent on the mean annual continental runoff, combined with a dependence on the mean annual air temperature through Arrhenius formalism (Walker et al. 1981). Continental weathering fluxes are spatially resolved at the FOAM general circulation model resolution (7.5° longitude by 4.5°l atitude). The calculated weathering fluxes also depend on continental bedrock lithology (granite vs basalt; Dessert et al., 2001;Oliva et al., 2003). Given the rapid dissolution rate of carbonate rocks, we assume that continental runoff water is at equilibrium with respect to calcite for each continental grid cell under the calculated CO 2 level. Sulphur dioxide reactions are computed in the atmospheric box using the parametric laws developed by Pierazzo et al. (2003) and Miles et al. (2004) . The 6 115 120 125 130 incoming radiative flux reduction caused by atmospheric aerosols is directly expressed by ΔT, depending on the sulphate aerosols loads (H 2 SO 4liq ) and sea ice thickness (for details, see Mussard et al., 2014). Recent simulations performed by Schmidt et al. (2016) are used to adjust the climate response of the model (see fig.2

electronic supplementary material).
One-third of the atmospheric aerosols are assumed to fall over continents, and the remainder go directly to the surface ocean. Over continents, sulphuric acid aerosols instantaneously react with the available carbonates and alter the weathering process (Donnadieu et al., 2006;Guidry and Mackenzie, 2000). The riverine input flux of sulphate ions into the oceans is calculated as the product of the runoff and ion concentrations. For oceans, H 2 SO 4liq primarily acts on the alkalinity budget and associated processes (e.g. the carbonate saturation state). Sulphate anions are also transported into the global ocean by ocean mixing within a few thousand years. In the absence of fully coupled simulations, including active oceanic dynamics, oceanic transport processes occurring in the biogeochemical model are represented by exchange fluxes between the ten-box model. These fluxes have been calibrated on present-day conditions and have held constant in all the simulations, given the lack of information about the changes in ocean mixing in the geological past.
To assess the oceanic environmental disturbance, long-term simulations (>100,000 yrs) are required. Since the length of simulations of general circulation models cannot exceed a few thousand years, the coupling procedure between FOAM and COMBINE is done through interpolations of the climate runs. Moreover, FOAM is used to generate an off-line catalogue of continental air temperature (T air ) and continental runoff (R) values for a wide range of pCO 2 values.
At each time step of the GEOCLIM calculation, and for each corresponding atmospheric pCO 2 , T air and R above each continental pixel are calculated through a linear interpolation procedure from the climatic catalogue.
In the GEOCLIM model, the oceanic biosphere is modelled as an organic carbon reservoir in the photic zone in each oceanic basin. Up to now, these reservoirs are fed by a flux of primary productivity, dependent on the input of phosphorus (the only nutrient modelled) inside the surface 7 140 145 150 155 ocean by upwelling or continental discharge. Biomass is then lost by decay and transferred to the deep ocean. In the present study, we introduce a highly simplified model for the biodiversity inside the photic zone. A fixed number of planktonic functional types (PFTs) n is prescribed (n = 50 in this study). These PFTs are represented by a mass of carbon (mol of C). Each PFT has its own physiological parameters and can interact with the others (primarily competing for resources ;Quere et al., 2005). The growth rate of each PFT, α i, is a function of the calculated surface temperature and pH according to the following equation: where br i is the birth rate, and α i is the maximal growth rate of species i, modulated by two Gaussian functions. The first one is centred on the optimum temperature of growth b 1 i with a tolerance of c 1 i , and the second is centred on the optimum pH of growth b 2 i with a tolerance of c 2 i . Additionally, the GEOCLIM model calculates the surface temperature T and the pH at each time step for each ocean basin. The optimal temperature of growth for each PFT, b 1 i , is fixed at the averaged value of the surface ocean reservoir for the pre-perturbation, with a random noise of ± 5°C. The same procedure is applied to the optimum pH, b 2 i , with the random noise fixed at ± 0.1 pH. The tolerances (standard deviation c 1 i and c 2 i ) are randomly set at a value within the interval [-4°C, 4°C] and [-0.3, 0.3] around the optimum temperature and pH.
The factor I i represents the potential competition for resources among the primary producers. For each species i, the existence of a competitive relationship with all the other species j is randomly fixed (yes or no). If there is no competition, I i is set to 1. Otherwise, it is calculated as the average of all competitive factors I i,j . These factors are calculated as a Michaelis-Menten function of the ratio of the biomass of PFT i and j. The death rate of each species is assumed to be proportional to their global biomass, with the proportionality constant being fixed randomly prior to the simulations. At this stage, we do not account for a trophic chain, including predators, nor t he 8 165 170 175 180 185 dynamics of ecosystems (the processes of evolution). All PFTs are thus primary producers. If the present version of the model does not explicitly simulate the biomass and behaviour of grazers, the effect of grazers is lumped into the death term of each producer, F p in , the nutrients flux (phosphorus).
To consider the role of carbonate ion concentration in calcification by marine organisms we used the CaCO 3 saturation state (Ω CaCO3 ). This assumption is based on positive correlation between Ω CaCO3 and the calcification rate, This link is defined as follows (Zeebe and Wolf-Gladrow, 2001):  (Stanlet and Hardie, 1999), Ω CaCO3 is largely determined by [CO 3 2-] . therefore the seawater pH. Using this property, if the ocean is still saturated (Ω>1) the calcification by marine organisms depends on the saturation state for calcite (see next section for details) or becomes null if (Ω<1). In this latest case the calcifying organisms are supposed extinct (their biomass is fixed to 0).

Boundary Conditions: The Late Maastrichtian
To quantify the K-Pg oceanic disturbance, the latest Cretaceous boundary conditions are applied to ECO-GEOCLIM. Climatic outputs are extracted from Maastrichtian simulations performed by Donnadieu et al. (2009). The land-ocean configuration (~68 Ma) is based on paleomagnetic data, hotspot tracks, geological constraints (Besse and Courtillot, 2002;Dercourt et al., 1993), the circular orbit of the Earth around the sun (eccentricity = 0), and the Earth's obliquity at 23.5° (this setting leads to an equal annual insolation for both hemispheres). In addition, the solar radiation is reduced by 0.6% (i.e. 1357 W/m²), according to the stellar evolution models (Gough, 1981). The tested atmospheric CO 2 values range from 160 to 2000 ppmv to cover all plausible atmospheric CO 2 scenarios for the K-Pg boundary (Caldeira and Rampino, 1990;Dessert et al., 2001). The presumed extent and spatial distribution of the Deccan Traps at the K-Pg boundary are based on the work by 9 195 200 205 210 Lefebvre et al. (2013). Plant cover is extracted from Donnadieu et al. (2009), where GEOCLIM was coupled to the dynamic vegetation model (LPJ; Lund-Postdam-Jena; Sitch et al., 2003). This coupling ensures constant coherence between climate and vegetation cover, which influences the silicate weathering but only through the runoff and temperature variables. The model does not account for other effects of vegetation on continental weathering (such as enhanced weathering by organic acids or the below ground CO 2 increase). Since FOAM does not consider the atmospheric chemistry of sulphur, land plants are neither affected by acid rains nor short-term cooling.
The Maastrichtian is considered a period of high carbonate production and deposition (Hay, 2004) with extant carbonate platforms and important chalk deposition (Kiessling et al., 2000(Kiessling et al., , 2003. In light of available data on change during the Phanerozoic, the real extent of tropical carbonate indicates that the latest Maastrichtian carbonate platforms (<30° latitudes) occupy ~16×10 6 km² (i.e., the modern areas of shelfal carbonate accumulation are ~0.6×10 6 km²; Walker et al., 2002).
Consequently, we updated the computation of the carbonate accumulation in shallow water (freef flux) as follows: freef = akcr * ocean surface * reefS* (Ω-1) 1,7 , where akcr represents the most dominant mode of carbonate production for shallow water (low/high value describes a Cretan/Neritan ocean, respectively), with reefS (areal extent of the tropical carbonate over the Phanerozoic period) as the constant 16×10 6 km². Thevolume of the coastal ocean is held constant at 4.86×10 6 km 3 .
In the neritic environment, carbonate production remains globally dominated by skeletal debris, such as large benthic foraminifera, rudists, bivalves, echinoderms, ectoprota, and algae with a small contribution from peloids and ooids (Kiessling et al., 2003;Steuber, 2002), while the open ocean was dominated by planktic foraminifera and calcareous nannoplankton. The ecological success of planktic calcifiers during the Mesozoic (Martin, 1995) induced a transition from a Neritan ocean, where the CaCO 3 production is driven by neritic organisms living on shelves, to a Cretan ocean dominated by biogenic pelagic CaCO 3 precipitation (Zeebe and Westbroek, 2003). To represent Cretan and Neritan oceans, we changed the percentage of the primary producers able to calcify in the photic zone. In the case of a Neritan ocean, 33% of the total biomass of the primary producers living on shelves are calcifiers (which corresponds to 17 PFTs as calcifiers and 33 as non-calcifiers; see fig

The K-Pg Events and Experimental Set-up
Using the Late Maastrichtian boundary conditions, the weathering of silicates and the simultaneous uptake of atmospheric CO 2 leads to stabilising the load of atmospheric carbon at 230 ppmv. This value is the result of elevated weathering rates due to emplacement of a large igneous province combined with the paleogeographic impact. Compared to usual estimates for the Late Cretaceous, 500 ppm < pCO 2 < 1500 ppm, our pre-crisis load of CO 2 appears too low (Royer, 2006). However the last CO 2 compilation shows that the atmospheric CO 2 at 65 Ma was close to 280 ppm, and may even be below that level (Foster et al., 2017) Moreover, a tropical climate prevails despite this apparently low load of carbon dioxide.
Simulated SSTs present an averaged value [50°S-50°N] of 18.6°C, accompanied by a wide latitudinal gradient, <15°C at mid-latitudes to >25°C close to the equator. Since this starting state provided us a description of the Late Cretaceous environment not so far for estimates proposed by proxies (Woelders et al. 2017), we preferred to leave the pCO 2 pre-crisis value unchanged, rather than tune the carbon cycle parameterisation (enhanced degassing and weaker weathering).
By fixing the latest Maastrichtian areas of shelfal carbonate accumulation at 16×10 6 km² (Walker et al., 2002), the alkalinity was depleted compared to present-day conditions; thus, the ocean was less saturated with calcite (Ω calcite ~ 2.2, instead of 4.5 < Ω calcite <5 for the modern ocean). the lysocline tends to be much shallower than present-day (2,7 instead of 4km) but no region of the surface ocean becomes undersaturated for calcite (Ω min polar regions = 2). This pre-crisis steady state implies that the latest Cretaceous ocean should be more sensitive to any process increasing the dissolved inorganic carbon without adding alkalinity.
To better explore the effect of coupling between environmental disturbances and biodiversity, we tested the effects of the Deccan Traps and/or bolide impact. To assess the effects of biological feedback in response to the Deccan Traps or bolide impact (or a combination of the two), we compared simulations where the primary production does not depend on environmental perturbations ('without biological feedback' experiments, Table 1) with simulations where the physiological parameters of species alter their efficiency in producing organic carbon ('with biological feedback' experiments, Table 1). In the absence of biological feedback, the amount of organic carbon produced is solely controlled by the available nutrients. Most of our runs used the Cretan ocean as the standard condition to fit palaeontological records; however, to highlight the role of pelagic calcifiers as primary producers, we also conducted experiments assuming a Neritan ocean. Table 1 summarises all runs performed in this study.

Eruptive Sequence of the Deccan Traps
The Deccan Traps outcrop over vast areas of west-central India and are currently formed by a 3500 m-thick lava pile. The eruptive sequence established in the Western Ghats, the thickest remnant of the main Deccan province, includes 30 major eruptive events and 41 single lava flows corresponding to a volume of ~253,000 km 3 (Chenet et al., 2009). After scaling to the total volume of the Deccan Traps (~2×10 6 km 3 ; Chenet et al., 2009), more than 75% of the lava volume is not considered, leading to an underestimation of the number of individual lava flows. Assuming a simple geometry for lava flows and accounting for the estimations of SO 2 and CO 2 concentrations in basalts (Self, 2006), Chenet et al. (2009) proposed that 15,000 to 35,000 Gt of CO 2 and 6800 to 17,000 Gt of SO 2 were released into the atmosphere, primarily during the emplacement of the main 12 275 280 285 290 295 province. We assume that sulphur species have reached the stratosphere due to the pulse-like degassing of the Deccan Traps (Kaminski et al., 2011), as suggested by akaganéite minerals (Font et al., 2017).
The emplacement of the Deccan continental flood basalts started with the eruption of the Northern Deccan province (Malwa plateau and Mandla area), a minor event in terms of volume, dated at 67.3 ± 0.3 Ma, lasting about 1.7 Ma (Schöbel et al., 2014). The main province started erupting during the C30n/C29r transition, continued during chron C29r, and finally, ended during chron C29n, straddling the K-Pg boundary (Chenet et al., 2007(Chenet et al., , 2008. The main province clearly represents the largest volcanic event during which at least 80% of Deccan lavas erupted. Chenet et al. (2009) suggested that the main erupting phase (Phase 2 according to Chenet et al. 2009) would have been as short as 200 kyrs. Recently, using U-Pb ages on zircon from lava flows collected in the Western Ghats (the main province), Schoene et al. (2015) concluded that the duration of emplacement lasted 0.7 Ma. Unfortunately, the absence of U-Pb ages out of the lowest and highest volcanic units did not allow a conclusion regarding the tempo of emplacement, notably the presence of a hiatus within the eruptive sequence. However, according to Schoene et al. 2019, this considered period covers the emplacement Polardpur and Ambenali formations (ie 66.1 to 65.9Ma), two among the most active phases of the Deccan traps. In addition, these two volcanic events encompass the K-Pg boundary (66.04Ma).
By deciding to focus on a shorter period of time of maximum activity (~200 kyrs) instead of studying 700 kyrs of eruption (Schoene et al. 2016, Sprain et al. 2019, we considered that the most important factor affecting the Earth's environment was the timing of eruption not the total duration of the emplacement of the Deccan traps. Indeed the emplacement of the main province (as all large igneous provinces) consists of individual lava flows and volcanic pulses made of successive flows erupting over a period that is too short to record the magnetic secular variation (< 100 yrs; Chenet et al., 2008Chenet et al., , 2009Courtillot and Fluteau, 2014). This result implies that the volcanic activity is interrupted by periods of quiescence, occasionally evidenced by thick weathered 13 300 305 310 315 levels (red bole). Following this evidence, we developed a sequence based on a pulse like activity instead of considering a linear degassing, a scenario appearing much more realistic to describe the large igneous province volcanism.
Based on constraints on time, volume, and fluxes of volatiles, we developed a stochastic model to obtain a pulse-like degassing scenario of CO 2 and SO 2 . We fixed the total volume of erupted magma at 2×10 6 km 3 , the total emission of CO 2 at 28,000 Gt, and the emission of SO 2 at 6,800 Gt. We assumed that the volume of lava flows and volcanic pulses follow a Gaussian distribution (m = 0, σ= 5000 km 3 ), meaning that only a few lava flows may exceed 10,000 km 3 . We

Gases Released and Solid Particles Ejected by the Chicxulub Meteorite
The Chicxulub crater (170 km in diameter) in the Yucatán peninsula is due to a collision with an asteroid of about 10 km in diameter. The sea floor at the target site is composed of a 3 km-thick sedimentary sequence of carbonate rocks (~35% dolomite and 30% limestone) and anhydrite (~30%) (Ward et al., 1995). To estimate the amount of CO 2 released into the atmosphere by vaporisation, we used the approach developed by Charnay et al. (2017). We tested three values of bulk density: 3320, 2490, and 1660 kg·m -3 to cover the entire range of cases (pure dunite with 0%, 25%, and 50% porosity, respectively). Assuming a decarbonation pressure fixed at 60 GPa (Martinez et al., 1995) and a mass fraction of CO 2 of 1% in asteroids (Charnay et al., 2017), the mass of CO 2 released by the impactor ranged from 1955 GtCO 2 to 1260 GtCO 2 . In line with these computations, we considered a mean value of 1400 GtCO 2 , a value in the range of Pierazzo's 14 325 330 335 340 345 estimate (1998) of ~1471 Gt. Our computational tests performed with the same assumptions predicted higher loads of SO 2 (~180 GtS or 360 GtSO 2 ). However, to be consistent with prior studies, while investigating the Chicxulub impact (Brugger et al., 2017), we prefer to use 100 GtS (or 200 GtSO 2 ) to represent the load of S released by the platform vaporisation. It is notable that additional sources of carbon due to burned biomass and remineralisation of the dying terrestrial vegetation were not considered.
The mass of ejecta produced by the impact, primarily silicates, is calculated using the approach proposed by Collins et al. (2005) and developed by Charnay et al. (2017). Assuming a density of 2490 kg·m -3 for the impactor, a total mass of ejecta distributed worldwide reaches 4.42×10 16 kg (2.95×10 16 kg with a density of 1660 kg·m -3 , and 5.86×10 16 kg with a density of 3320 kg·m -3 ). The weathering of ejecta, a net sink of CO 2 , is calculated following Sleep and Zahnle (2001). In the absence of modules devoted to include following complex processes nitrogen oxide formation due to atmospheric shock waves (Tyrrell et al., 2015), ozone depletion, and bedrock contamination (Ganino and Arndt, 2009), for example, were not considered. . Since the residence time of water in the atmosphere, it would not have exceeded a few months; thus, we do not consider its effect further.

Time Coincidence Between Deccan Traps and the Chicxulub Meteorite
Recent high-precision 40 Ar-39 Ar dating of the K-Pg boundary in Hell Creek (Montana, USA) and of the Chicxulub (Yucatán, Mexico) impact ejecta shows that these two events are time-coincident within ~32,000 yrs precision at ca. 66.04 Ma (Renne et al., 2013). The discovery of the oldest spherules in the vicinity of the Chicxulub astrobleme in layers is attributed to C29r (between C30n and C29n), suggesting that the impact occurred at least 10 kyrs later than the main phase of the Deccan Traps (Keller, 2014).
Nevertheless, the sudden outburst of Deccan eruptions, accounting for more than 70% of the main-phase eruption, occurred within ~100,000 yrs or less of the K-Pg boundary (Renne et al.,15 355 360 365 370 2015), suggesting a longer interval. Without precise time constraints, to enhance asteroid effects on the environment, we assumed that the Chicxulub impact struck the Earth's surface 100 kyrs later than the inception of the Deccan main province volcanism. It is noteworthy that a sensitivity experiment with a different time coincidence (10 kyrs instead of 100 kyrs) presents similar effects (ΔSST~0.1°C, see electronic material fig. 5).

Chicxulub Impact Modelling
Over the five years following the impact, the modelled SST drops by 8°C in the low latitudes despite the rise in CO 2 (ΔpCO 2 = +120 ppmv, Fig. 1) caused by rock vaporisation (especially carbonates). In response to sulphate aerosols (Robock et al., 2009;Schmidt et al., 2016), the cooling trend continues for at least three decades after the impact. Because the S-aerosol particles are progressively washed out from the atmosphere, the cooling effect is progressively outweighed by the warming provided by the pCO 2 due to its longer residence, which causes post-impact warming of ~1.5°C after a century (Fig. 1). Due to the ice-albedo feedback, the warming of polar SST is delayed by about 100 yrs and does not exceed +0.5°C.
The post-impact CO 2 evolution is highly sensitive to biological feedback. Succeeding a brief increase in CO 2 , pCO 2 decreases within 1 kyr in the absence of feedback (A-Cretan-2), while several 10 kyrs are needed when feedback is considered (A-Cretan-1). As long as S-aerosol amounts remain high, a transient shutdown of the primary production is simulated when biological feedback is considered (Fig. 2a). The primary productivity briefly increases within the first years after the end of cooling, a phenomenon that could be related to the phosphorus accumulation from previous years. This transient peak can be considered a 'swan song'. After the removal of the excess phosphorus, the primary productivity decreases and reaches its lowest value around a decade after the impact (Fig. 2a). Between 50 and 200 years after the impact, the consumption of nutrients is less 16 380 385 390 395 efficient due to environmental changes ( Figs. 1 and 3), and the phosphorus content increases (Fig.   2b).
Driven by oceanic disturbances, despite the presence of nutrients (Fig. 2b), primary productivity and oceanic biomass both decline (Fig. 4), especially the calcifying primary producers, due to acidification of the surface ocean (Fig. 3). The recovery of the primary productivity happens over a period of a few centuries (Fig. 2) when pre-impact conditions tend to be restored. In response to the low percentage of residual living biomass (Fig. 4) in the aftermath of the impact, the burial of organic carbon is weakened; thus, only weathering processes efficiently reduce pCO 2 .
Consequently, by including biological species' inherent properties, environmental perturbations lead to creating a sluggish organic carbon cycle characterised by a delayed response of reducing pCO 2 .
However, whatever seawater conditions are tested (Cretan or Neritan), the surface ocean remains oversaturated with calcite (Ω calcite ~ 1.3; Fig. 3c). A puzzling behaviour is observed a few decades after the asteroid impact. Despite the ocean acidification, the model predicts an increase in carbonate precipitation (Fig. 3d), while the surface ocean remains depleted of CO 3 2-(1.3 < Ω calcite < ~1.5; Fig. 3c), a result that is contradictory to the theoretical view, assuming a strict dependence on carbonate deposition with Ω calcite . In fact, because carbonates are bio-induced (see Section 2.2 for details), marine calcifier organisms benefit because all primary producers of abundant nutrients remain available in the surface ocean. Consequently, during the first century after impact, carbonate deposition is dominated by primary productivity rather than by Ω calcite . With For a better understanding of how oceanic disturbances lead to extinctions among calcifying and noncalcifying species, we explored variations of living marine organisms in a Cretan ocean before, during, and after the Chicxulub impact. We developed the following relationships for each PFT in the photic zone, where primary producers are abundant (surface boxes). The present relationship is based on specimen abundance (SA) for each PFT (full details given in footnote 1 ) and provides an apparent extinction pattern (i.e., apparent species extinction or aSE).
Based on our assemblage of marine PFTs, A-Cretan-1 case illustrates that, if the overall disturbance is sufficient to trigger a biomass reduction of 60% for calcifying organisms (40% for noncalcifying; Fig. 4), this results from the collapse of a few dominant calcifying PFTs that were providing the most biomass. Furthermore, biotic losses for primary producers would not have exceeded 20%, assuming an SAT close to 0.1 (i.e., the abundance of living specimens decreases to 10% of the pre-crisis value). Despite these drastic perturbations, this result suggests that no major extinction occurs when we combine an instantaneous release of 1400 GtCO 2 with 200 Gt of SO 2 .
Additional sources of carbon due to burned biomass and remineralisation of the dying terrestrial vegetation may enhance this perturbation however but SO 2 emissions cannot rise above 250Gtthreshold of the reversible cooling (see section 4.1 to see sensitivity experiments).

1
We defined the mass extinction ratio based on specimen abundance. We supposed that each type of organism remains characterised by a body shape and mass held constant through time. Biomass (t = 0) and (t) correspond to the initial and computed biomass at time t, respectively. The specimen abundance threshold (SAT) represents an arbitrary value to describe the apparent presence of specimens in the fossil record. In this study, the arbitrary SAT ranges from 0.2 to 0.01. This range is used to represent the complex set of factors distorting the observed extinction record for each species (the Signor-Lipps effect). According to this assumption, when SAT is fixed to 0.01, the extinction occurs when the apparent abundance of living specimens reduces to 1% of the pre-crisis value. Based on these equations, species can reappear later in the geological record if the biomass of the species rises. 18

Deccan Traps Modelling
To simulate the effects of the Deccan Traps, we computed a sequence of 200 kyrs with pulses following the stochastic distribution as described in Subsection 2.3.1. The response to a single volcanic pulse presents tendencies similar to those simulated for the Chicxulub impact but with a weaker amplitude. Both cases show an abrupt cooling over several decades, driven by S-aerosols, followed by warming caused by enhanced pCO 2 after aerosol scavenging. However, through time, the load of CO 2 accumulated by degassing by the Deccan Traps overcomes the rise in CO 2 triggered by a 10-km impactor, +250 < ΔpCO 2 < +470 ppmv (instead of +120 ppmv for the Chicxulub impactor; Figs. 5 and fig.1 in electronic supplementary material).
The primary driver controlling the rise of pCO 2 is the degassing itself; however, biological feedback suggests an external mechanism able to amplify the initial perturbation. The mean response obtained with the D-Cretan-1 scenario presents a rise of SST by 3.5°C (~23.5°C) between 50°S and 50°N, whereas the warming does not exceed 2.5°C in the absence of biological feedback (D-Cretan-2). This finding demonstrates that the same eruption scenario generated variable environmental disturbances, depending on the properties of the oceanic ecosystem (Fig. 5).
Nevertheless, the style of an eruption has a significant influence as well. Examining the linear degassing case (Fig. 5), volcanic sulphur products are washed without disturbing the long-term warming provided by the CO 2 (the warming is just slightly weaker due to the continuous S-aerosol formation; D-Cretan-linear) whereas the opposite behaviour is found when the eruptive sequence includes a succession of pulses (the long-term warming is interrupted by a succession of cooling events; D-Cretan-1). The slowing of the weathering process caused by the aerosol-induced cooling is highly dependent on the assumed quiescence period between two pulses. For instance, using a quiescence period twice longer D-Cretan-long-degaz), the excess CO 2 is removed from the atmosphere because silicate weathering becomes efficient enough to limit long-term warming (i.e. 19 455 460 465 470 self-regulation of the Earth's temperature; Walker et al., 1981). Consequently, if the eruptive sequence is intense enough, the long-term warming should be briefly interrupted by cooling events lasting 10 to 100 years. Otherwise, the Earth's climate should be characterised by limited long-term warming due to the efficiency of silicate weathering. Despite the importance of emitted SO 2 (D-Cretan-noSO 2 ) and aerosol formation in the Earth's climate by limiting silicate weathering (Mussard et al., 2014), the intense build-up of CO 2 caused by biological feedback illustrates that organic carbon more efficiently destabilises the global carbon cycle than inorganic processes over the time scale of a thousand years.
Over the first 100,000 years of volcanic activity, what happens at the surface vastly affects the deep ocean, and both parts present the same tendency, although the shift in deep water towards lower pH is greater than the surface ocean ( Fig. 6a). At the opposite end, over the next 100,000 years, the acidification of the surface ocean stops (Fig. 6b), while the pH of the deep ocean still decreases (Fig. 6a). Acidic seawater affects the deep and surface oceans in varying ways. Most of the calcite is dissolved at depth (Fig. 6c), forcing the lysocline to rise close to 1500 m (Fig. 6d). The The results from oceanic disturbances modelled in D-Cretan-1 run is a biomass reduction of about 90% for calcifying organisms and 65% for noncalcifying organisms (Fig. 7). Conversely to the asteroid case (A-Cretan-1), the biomass fall simulated with Deccan Trap volcanism (D-Cretan-1) corresponds to a global collapse, whatever the considered PFTs. Calcifying PFTs are more affected because their primary productivity efficiency is affected by the cumulative effect of abiotic factors varying (ΔSST, ΔpH, and calcite saturation), whereas no calcifiers remain only affected by SST and pH changes. It is noteworthy that most parts of the biomass lost occurred on a time scale ranging from 10,000 to 100,000 years, suggesting that an eruptive sequence that will last longer is not necessarily associated with more drastic environmental changes and biomass crises. If a biomass disturbance seems sufficient to qualify this event as a massive collapse, the patterns of extinction for calcifying and/or noncalcifying PFTs, both ranged from 50% to 80%, suggest a mass extinction event in terms of marine biodiversity as well.
To specifically investigate the effects of seawater pH on calcifying and noncalcifying producers, we conducted an additional run in which pH fluctuations have no effect on the marine biosphere (D-Cretan-3). Compared to standard conditions (D-Cretan-1), after 100,000 years of volcanic activity, the phosphorus content becomes twice lower, while the primary productivity is only reduced by 40% instead of 60% (Fig. 8). However, in two cases, we observed a clear inflexion over the next 100,000 years with the re-increase of primary productivity. This tipping point is crossed when the ΔpH exceeds 0.1. Since most of the marine PFTs appear tolerant to Δ pH of ~0.1 ( Fig. 3 SI), the organic carbon cycle will not be destabilised below this threshold of fluctuation. As long as changes in seawater pH do not exceed 0.1, the primary production remain nearly unaffected, and the productivity curves for both experiments look similar through the first 60-70 kyrs.
After this timespan, we observe that the primary production and phosphorus tendencies appear less pronounced for the D-Cretan-3 experiment (only driven by variations of SST). Finally, in response to a sturdier rate of organic carbon burial linked to a more abundant biomass in oceans, 21 510 515 520 525 the D-Cretan-3 run presents a lower pCO 2 (ΔpCO 2 = +420 ppmv) compared to standard conditions (ΔpCO 2 = +470 ppmv) at the end of the Deccan Traps.

The K-Pg Mass Extinction Modelling
We now combine the 200 kyr-long eruptive sequence of the Deccan Traps and the Chicxulub impact, which collided with the Earth 100,000 years after the onset of the degassing of Poladpur In addition to climate changes, atmosphere-ocean exchanges lead to a pH decrease in the global ocean by 0.2 for the surface ocean (from 8.0 to 7.8) while the deep ocean pH drops by 0.4 (from 7.8 to 7.4). Without requiring any additional perturbations, these environmental disturbances lead to a dramatic decrease in calcifier organisms (~90% of their biomass and ~80% of their apparent biodiversity, assuming an SAT fixed at 0.1), while noncalcifying organisms appear less affected (biomass and biodiversity are reduced by about 65% and 60%, respectively) . It is noteworthy that mass extinction patterns (Fig. 9e) are affected by the Chicxulub impactor, and the mass extinction ratio of calcifying and noncalcifying species only differs after the impact, suggesting that very rapid marine acidification can be considered an effective trigger for explaining the severity of extinctions among pelagic calcifying species compared to non-calcareous taxa.

Discussion
Here, we compare our modelling results with prior studies and geological records to validate the mechanisms described in Section 3.

Chicxulub Impactor
According to Brugger et al. (2017), the mid-tropical SST underwent a similar cooling (ΔSST = -8°C in this study compared to -6°C), whatever the stratospheric residence time assumed for S-aerosols.
Both studies used the same amounts of released gases, but the S-aerosol computation is governed by different processes. In ECO-GEOCLIM, the S-aerosol loads are derived from chemical reactions depending the gas amounts with the effect on the climate predicted using a parameterisation based on that by Schmidt et al. (2016). Conversely, CLIMBER3α follows the assumptions by Pierazzo et al. (2003), but the S-aerosol load is simulated by a reduced solar flux. However, in the absence of direct measurements (high loads of S-aerosols not yet observed), considering the complex interplays of sulphur species in the atmosphere (aerosol size; Timmreck, 2012), the scavenging issue (Ohno et al., 2014), and the altitude of ejection (Kaminski et al., 2011), both estimates of climate change can be considered consistent with the existing data. In addition to these processes, the cooling could be more marked using a mixed-layer oceanic model for computing SST (this study), with the heat transport computation being less accurate than in a fully dynamic ocean model (Brugger et al., 2017).
The reconstructed SST based on TEX 8 6 palaeothermometry for sediments deposited within 100 years following the impact (Vellekoop et al., 2014) also shows a sharp drop in tropical coastal seawater temperature (ΔSST max = -7°C). The same section (Brazos River, Texas) also reveals a prolonged warming following the 'impact winter', where the post-impact SST is 1°C to 2°C warmer than the pre-impact values. Compared to our model results for the Chicxulub impact, the best match is observed with the A-Cretan-1 run. Indeed, the data and model results show the same 23 560 565 570 575 580 sequence of events: intense cooling followed by moderate warming. Without excluding some issues in the recording of such short-lived events (reworked materials), palaeothermometry measurements on the Brazos River section also fit well with the AD-cretan-1 run (Table 1), where distinct cooling phases could be potentially associated with volcanic pulses in addition to the Chicxulub impact.
Whatever the scenario considered, our model failed to reproduce massive acidification ( Fig. 3) that would be able to create surface waters that are undersaturated with calcite, as defended by Ohno et al. (2014). If the S-aerosol scavenging within a few days may be invoked to solve this issue, this assumption also implies small amounts of stratospheric S-aerosols and a lifespan of Saerosols that is too short to explain the Earth's cooling (Vellekoop et al., 2014).
Following this conclusion, we explored an alternative solution by carrying out a set of experiments where we test the Earth's system response to massive injections of sulphur (D'Hondt et al., 1994;Ohno et al., 2014;Tyrrell et al., 2015). Holding the mass of CO 2 released by the impact constant, sensitivity tests revealed that surface waters become very briefly (< 1 year) undersaturated for calcite (Ω calcite ~ 0.95), with a modest increase of injected sulphur (175 GtS instead of 100 GtS).
Very massive injections can help to extend the length of the surface ocean undersaturation of calcite, but such high fluxes lead to surface cooling of the Earth. For example, the tropical SST [50°S-50°N] declines to 0°C for 1 year in response to an injection of 250 GtS (surface ocean pH = 7.6, Ω calcite~0 .75) or 10 years with 770 GtS (surface ocean pH = 7.4, Ω calcite < 0.5). Beyond this amount of sulphur, an irreversible cooling occurs due to sea ice spreading and snow extending.
Within a few years, the Earth becomes entirely covered in ice (snowball Earth state). Although climate instability thresholds due to ice-albedo feedback may be model-dependent, the first-order conclusion rejects the existence of very massive injections of sulphur, a conclusion that is in agreement with findings by Tyrrell et al. (2015). Consequently, based on AD-Cretan runs (Fig. 9), we defend that the likely solution will be a combination of long-term acidification provided by CO 2 24 590 595 600 605 (Deccan Traps) coupled with short-term acidification resulting from the rainout of sulphate aerosols (Chicxulub).

Deccan Traps Emplacement
Using δ 18 O measured in foraminifera, the Late Maastrichtian (chron C29r -planktic foraminiferal zones CF2, then CF1 -last biozones before the K-Pg boundary) was originally identified as a tropical climate (Li and Keller, 1998). This view was challenged by recent paleoclimate estimates that suggest a succession of warm-cool events rather than a single global warming event (Punekar et al., 2014). These new data lead to a reassessment of the K-Pg boundary, which may be defined as follows (Keller, 2014): the CF2 zone (spanning ~120 kyrs) characterised by an intense warming close to 3°C to 4°C (Li and Keller, 1998), followed by the CF1 zone (~160 kyrs), a warm fluctuating period because the warming trend ends in the middle of the CF1 zone, and the Earth's temperature starts to decline. Despite this complexity, all experiments simulating Deccan Traps emplacement (with or without Chicxulub impact) gather all explain some of the major features presented in the geological record: 1) the combination of warm-cool events, 2) warming duration (150-200 kyrs) fitting with foraminifera biozones, and 3) the SST rise in close agreement with proxies ~3.5°C, which is roughly equivalent to 7°C in high latitudes, according to Tobin et al. (2017). The model-proxy agreement emphasises that the emplacement of Deccan Traps seems to be the main driver causing climate changes over the K-Pg boundary. Tobin et al. (2017) demonstrated the relevance of pulse-like degassing to explain climatic perturbation and revealed a striking correlation between warming and emitted CO 2 . In their attempt to reproduce 3°C to 6°C of global warming (7°C for polar regions) deduced from oxygen isotopes (Barrera and Savin, 1999;Li and Keller, 1998), Tobin et al. postulated that 35,000 Gt CO 2 must be injected over 200 kyrs for rising pCO 2 by ~400 ppmv and that it warms the Earth's surface by ~3°C.
These estimates were obtained with a box model describing the biogeochemical cycles (GEOCYC, a model derived from GEOCARB) in which the load of atmospheric CO 2 is solely driven by the 25 615 620 625 630 enhancement of silicate weathering in response to Earth's heating, with feedback precluding the accumulation of CO 2 into the atmosphere. Considering the biological and SO 2 feedback, we revised this issue because both mechanisms act as a limiting factor of the carbon consumption. Our model simulates more intense global warming of 3.5°C and +470 ppmv for a lower gas release compared to that found by Tobin et al. (2017) (28,000 GtCO 2 and 6800 GtSO 2 ).

The K-Pg Calcareous Plankton Diversity Conundrum
The simulated biomass of calcifying and noncalcifying plankton (Fig. 9d) clearly captured the palaeontological pattern of a more important turnover of calcareous plankton (calcareous nannoplankton and planktic foraminifera) in comparison to siliceous and organic plankton (diatoms and dinoflagellates; e.g. Macleod et al., 1997). This difference in plankton type response to the environmental perturbation clearly lies in the carbonate system and pH degradation, which is more detrimental for calcifying organisms.
Unfortunately, quantitative and semi-quantitative studies on calcareous nannofossil accumulation rates show no decrease before the K-Pg boundary (Bernaola and Monechi, 2007;Henriksson, 1996;, suggesting that the calcareous nannoplankton productivity (i.e., organic production by calcareous nannoplankton approximated by the accumulation rate of calcareous nannofossil) was not affected by the environmental changes (pH, Ω calcite , and SST) induced by the onset of Deccan Traps Phase 2. On the other hand, the mass extinction event in relation with the Chicxulub impact clearly decreased the calcareous nannoplankton productivity (Bernaola and Monechi, 2007;Henriksson, 1996;Hull et al., 2011). This singular feature of nannofossils leads the authors to conclude that the Chicxulub impact has a more important influence than the Deccan Traps activity (e.g. Bernaola and Monechi, 2007;Bown, 2005;Gardin, 2002).
This apparent discrepancy between nannoplankton fossils and our model runs may be explainedby Fig.9e. According to these results, the Deccan traps should be responsible of the 26 640 645 650 655 collapse of primary producers while the abrupt short-lived acidification caused by the Chicxulub impact may explain the severity of the extinction among pelagic calcifying species.

The K-Pg Mass Extinction
With the K-Pg mass extinction, a collapse of the organic flux from the surface to the deep sea is described due to the decrease in plankton organic production (Strangelove Ocean Hypothesis; Zachos et al., 1989). More recently, a reduction of export of production to the seafloor by a change in organic matter recycling was proposed (Living Ocean Hypothesis; D' Hondt et al., 1998).
Nevertheless, these studies were challenged by recent findings. Indeed, the benthic foraminifera did not suffer significantly from the mass extinction event (Alegret et al., 2012;Kaiho, 1992), whereas the organic abundance of algal steranes and bacterial hopanes alongside the recovery of the δ 13 C org and δ 15 N org ratio would suggest a brief (less than a century) reduction in the algal primary productivity (Sepúlveda et al., 2009). A primary productivity proxy, the biogenic barium accumulation rates would suggest a more complex pattern: on one hand, a rapid recovery to no effect in the Central Pacific Ocean or neritic regions of the Atlantic Ocean and, on the other hand, a depressed organic flux in the Southern Ocean, Indian Ocean, and northeast and southwest Atlantic Ocean (Hull et al., 2011).
Our simulations propose some plausible explanations for this regional variability. nutrients to the surface may remain relatively unaffected, with the load of nutrients overcoming the low efficiency to synthesise organic matter.
Among the biologic oceanic community, the mass extinction event has not affected all the groups the same way (Macleod et al., 1997). The coccolithophores and associated incertae sedis, altogether defined as calcareous nannoplankton, were particularly affected by the mass extinction event with a loss of 93% of the diversity (Bown, 2005), corresponding to the most important diversity loss of this group since its first occurrence in the Late Triassic (Bown et al., 2004). The pattern of this turnover is particularly interesting. Although the diversity loss seems instantaneous on the geological scale, before the K-Pg boundary, the species richness and relative abundances fluctuate in relation to climatic changes at the end of the Maastrichtian. Most notably, the episode of warming, lasting ~200 to 250 kyrs before the K-Pg boundary, has induced an increase in relative abundance of warm-water species, extending the latitudinal range of presence of those species, and was marked by a slight decrease in species richness Thibault and Husson, 2016).
This global warming seems to be associated with the beginning of emplacement of Deccan Traps Phase 2 (Petersen et al., 2016;Thibault and Husson, 2016). Nevertheless, as noted by Thibault et al. (2015), those changes in species richness are relatively weak compared to the diversity and absolute decrease in abundance of calcareous nannofossils in the sediment at the K-Pg boundary (Bernaola and Monechi, 2007;Henriksson, 1996;Hull et al., 2011).
Immediately following the mass extinction, the calcareous nannoplankton community is still dominated by the Cretaceous survivor taxa, which is rapidly replaced by the newly evolving Paleogene species (e.g. Bown, 2005). The switch in the community is associated with a decrease in the size of the coccoliths. The broad pattern is the loss of large Cretaceous species, which are replaced by smaller Paleogene species (Herrmann and Thierstein, 2012). Along with the absolute abundance decrease, this exacerbates the decrease in the carbonate accumulation rate. Nevertheless, in detail, this pattern is challenged by the rapid occurrence in the Paleogene of large coccolith 28 690 695 700 705 710 (Bernaola and Monechi, 2007;Bown, 2005).The Cretaceous survivors are supposed to be rstrategist species (opportunistic), which are more frequently in high latitudes and neritic environments (Bown, 2005). Those species were more likely to survive a mass extinction event because of their capacity to live in unstable environments with changes in temperature, nutrient concentration, or even Ω calcite on a short time scale (days, seasons, or years).
Looking at all of these studies, if some aspects of our results are in accordance with the observed palaeontological pattern, our initial assumptions appear still too simple to capture the complexity of the patterns of K-Pg extinction. Consequently, we highlight in the following section the limitations of this study to present the future improvements.

Limitations
Our simulations provide an overview of possible patterns of marine primary productivity and biodiversity collapse in response to climate changes and oceanic acidification during a short window (200 kyrs) of the K-Pg transition (~1 M yrs). Our simulations suggest that emissions of SO 2 and CO 2 are first-order triggers of the K-Pg environmental crisis. Considering that short-term processes reveal that the initial perturbation (gas loads of CO 2 and SO 2 ) may be amplified by two major factors: 1) the type of eruption over the course of 200 kyrs and 2) the inherent biological parameters of marine organisms. Unfortunately, the time computation has led us to limit the number of factors tested. To advise the reader about implications of these choices, the cautions are briefly presented with respect to the limitations of this study.
To assess the potential effect of the Deccan Traps, we conservatively assume that 3.4 Tg of SO 2 and 14 Tg of CO 2 were degassed for every cubic kilometre of erupted lava. These values are held constant over 200 kyrs. With this approach, we omitted 1 Tg of chlorine per cubic kilometre yield by the Traps (Self et al., 2008). In addition, eruption parameters, such as the effusion rate and tempo, have been considered important to explain environmental effects (Rader et al., 2017).
Unfortunately, the lack of information about these parameters has led us to build the eruptive 29 720 725 730 735 sequence using a stochastic approach. Despite the development of a pulse-like emission, the absence of a realistic degassing scenario at the K-Pg boundary implies that major environmental features occurring before or after the Deccan Traps cannot be captured (i.e., pre-event warming accompanied by the rise of M. Murus; Thibault et al., 2017, the post-Traps overshoot or the absence of a progressive warming 100 kyrs before the K-Pg boundary).
In addition, we addressed sensitivity to climate change by computing inherent parameters with random functions (temperature and pH dependence), without considering a species distribution model or the ocean dynamics. Consequently, the presented results are model-dependent. The primary producer collapse is driven by the set of assessed species. To balance the lack of data to establish the tolerance patterns of species, only simulations including a high number of species (1000 or more) and/or different sets of marine biospheres would be more accurate to capture all patterns recorded in the fossil database. The K-Pg boundary is well known for extinctions of marine species occupying the upper levels of the food chain, such as ammonites, belemnites, or marine reptiles. This study does not include any predator-prey model. In the absence of such links, we cannot explore potential extinctions of consumers due to cascading effects in the food web. This aspect is beyond the scope of this study, but this feedback should be tested in future models to determine why calcareous plankton were more severely affected by the mass extinction event than other calcareous benthic organisms.
The simple model developed here to describe oceanic primary productivity is not yet an ecological model. However, for the first time in a long-term climate-carbon model, it introduces a retroaction between marine biodiversity and the long-term evolution of biogeochemical cycles. This study introduces the development of more complex trophic chains featuring several levels of consumers.
Regarding the Chixculub impact, we do not explicitly consider the change in the light influx (i.e impact winter assumption by blocking out of sunlight) because marine organisms may present tolerance to a decline of light flux by increasing the size of organic structure absorbing photons 30 745 750 755 760 (Geider et al. 1997). Consequently, the change of light influx has been considered by computing the primary productivity as a function of SST (see section 2.1). This assumption is based on measurements of photosynthetic rates derived from satellites (Behrenfeld and Falkowski, 1997validity range -1 to 29°C).
Finally, we cannot clearly assess the K-Pg transition without more precisely exploring processes occurring over continents. Traces of acidification in lacustrine sediments (Font et al., 2016) and strong weathering related to volcanic pulses (Gertsch et al., 2011) suggest that the neutralisation reaction is not verified everywhere (i.e., sulphuric acid reacting with limestone), as defended by Maruoka and Koeberl (2003) and supported by the very low extinction rate of freshwater vertebrates (Sheehan and Fastovsky, 1992). Consequently, very strong acidification may have occurred in rivers or restricted zones of the coastal ocean. In addition, evidence of very drastic climate changes over continents (Tobin et al., 2014) may also enhance extinction rates in the continental domain. Continental species should disappear more rapidly than marine organisms.

Conclusion
Using a coupled biodiversity-climate-carbon numerical model (ECO-GEOCLIM), we explored the triggering mechanisms of extinctions through the K-Pg boundary by quantifying the environmental effects on primary productivity, biomass, and biodiversity. We discovered two plausible pathways      fig.4 ). According to the AD-Cretan 1 experiment, the time coincidence between the traps of the impact of Chicxulub is fixed in half of the Deccan traps. The left/right panel presents SST before/after the asteroid fell in Yucatan Peninsula (black circle).
Corresponding oceanic disturbances are detailed in the fig.9