2.1 Sediment cores

Core SN^∘4 (27^∘59^′46.4^′′ N, 12^∘32^′40.6^′′ W) was retrieved in the Tarfaya Basin (southern Morocco), 40 km east of the town of Tarfaya, close to the road to Tan-Tan. The 350 m long marine sediment succession in Core SN^∘4, which provides an expanded record of OAE2, was deposited in an outer shelf setting within a subsiding basin on a passive margin during the Late Cretaceous. Late Albian to Turonian stable isotope and geochemical records were previously presented by Beil et al. (2018) and Scholz et al. (2019). High-resolution XRF scanner, bulk sediment stable isotope and carbon records across the onset of OAE2 were provided by Kuhnt et al. (2017). Here, we complement these records with phosphorus speciation data and time series analysis of logging natural gamma ray and high-resolution XRF scanner elemental distribution data across OAE2.

Cores LB1 (43^∘14^′37^′′ N, 5^∘34^′12^′′ E; 70 m long) and LB3 (43^∘14^′42^′′ N, 5^∘34^′52^′′ E; 56 m long) were drilled near Roquefort-La-Bédoule, 16 km southeast of Marseille in southern France. During the Aptian, the coring sites were located within an isolated intra-shelf basin, the South Provence Basin, on the northern Provençal carbonate platform, where an expanded succession of marine sediments including OAE1a accumulated. Details of the drilling operation were published by Flögel et al. (2010), lithologic descriptions, biostratigraphy and intermediate-resolution carbonate and stable isotope measurements were provided by Lorenzen et al. (2013) and Moullade et al. (2015). Here, we present new high-resolution bulk carbonate isotope data across the precursor, onset and peak phases of OAE1a as well as time series analysis of logging natural gamma ray and XRF scanner elemental distribution data from a spliced composite record of LB1 and LB3. The tie point between the two cores is at a depth of 49.4 m (Section 34, 3 cm section depth) in LB3 and 10.83 m (Section 8, 95 cm section depth) in LB1, below bed 170, identified by Moullade et al. (2015) as a prominent correlating feature. The tie point was selected at the base of a prominent, extended maximum in XRF-scanner-derived log(Terr∕Ca) (Sect. 2.3) identified in both cores.

2.2 NGR logging

Wireline borehole logging of the SN^∘4 drill hole was carried out by Geoatlas (Laayoune, Morocco) using a Century geophysical logging system with a natural gamma ray (NGR) sensor. Detailed information on NGR logging of SN^∘4 is given in Beil et al. (2018). For borehole logging of the LB1 and LB3 drill holes, an Antares Aladdin logging system with GR5 sensor probe was deployed. All logging operations were conducted shortly after completion of the drilling operation in boreholes without metal casing, and NGR data are presented as American Petroleum Institute radioactivity units in counts per second (cps) with vertical resolutions of 0.1 m for SN^∘4, 0.05 m for LB1 and 0.025 m for LB3. The intensity of natural gamma radiation is predominantly influenced by the concentration of three different elements: potassium (K), uranium (U) and thorium (Th). All three elements are bound to clay minerals, with uranium also adhesively enriched in organic matter. In environments with low terrigenous input and high organic matter accumulation as in Core SN^∘4, NGR is predominantly controlled by the concentration of organic matter. By contrast, clay accumulation mainly controlled NGR in Cores LB3 and LB1, where the organic matter content of the sediment is low.

2.3 XRF scanning and line scan imaging

Detailed descriptions of the XRF scanning of Core SN^∘4 are provided in Kuhnt et al. (2017) and Beil et al. (2018). Sections were scanned with a second-generation Avaatech X-ray fluorescence scanner at Kiel University. Surfaces were covered with a 4 µm thick Ultralene foil after cleaning with fine-grained sandpaper. Data for this study are from the 10 kV dataset scanned with a spatial resolution of 1 cm (vertical slit of 1 cm, horizontal slit of 1.2 cm) measured with 10 kV, 750 µA, no filter and 10 s acquisition time. Sections of Cores LB1 and LB3 were polished with fine-grained sandpaper and covered with 4 µm thick Ultralene foil prior to scanning with the Avaatech XRF scanner at Kiel University. The XRF data from the 10 kV dataset presented in this study were measured with a spatial resolution of 1 cm with the tube setting of 10 kV, 250 µA and no filter with 15 s acquisition time.

Raw XRF spectra from all three cores were converted with the iterative least-square package Win Axil batch (Canberra Eurisys) and the 10 kV Kiel model into element-specific area counts. The Jai CV-L 107 3 CCD colour line scan camera installed in the Avaatech XRF scanner at Kiel University was used to obtain core images with a downcore resolution of 143 ppcm for Cores SN^∘4, LB1 and LB3. We used the logarithmic ratios of elemental counts to eliminate XRF-scanning-specific effects such as the matrix effect, grain size effect or variations in rock density (Weltje and Tjallingii, 2008).

Log(Terr∕Ca) is used as a proxy for terrigenous, clastic material vs. marine biogenic carbonate. Elements of typically terrigenous origin (Terr) detected with the 10 kV settings are aluminium (Al), iron (Fe), potassium (K), silicon (Si) and titanium (Ti) (e.g. Peterson et al., 2000; Calvert and Pedersen, 2007; Mulitza et al., 2008; Tisserand et al., 2009; Govin et al., 2012). The variability of Fe closely matches the variability of the other terrigenous elements (Sect. S2), thus appearing not to be influenced by changing redox conditions. Calcium (Ca) is of marine origin, mainly from calcareous nanoplankton and microplankton. Log(K∕Al) is widely used to characterize the composition of clay mineral assemblages (Weaver, 1967, 1989; Niebuhr, 2005). The ratio is primarily controlled by variations in the amount of potassium-rich illite and thus reflects the intensity of physical and chemical weathering in the source area (Calvert and Pedersen, 2007).

2.4 Stable isotopes

Organic and bulk carbonate stable isotope data from Core SN^∘4 were compiled from the high-resolution record over the onset of OAE2 (Kuhnt et al., 2017) and from the lower-resolution record over the entire core (Beil et al., 2018). Bulk carbonate samples were analysed with the Finnigan MAT 251 and MAT 253 mass spectrometers at the Leibniz Laboratory for Radiometric Dating and Stable Isotope Research at Kiel University. The accuracy for carbon isotopes was ±0.05 ‰ and ±0.08 ‰ for oxygen isotopes. Organic carbon stable isotopes were analysed at the GeoZentrum Nordbayern (Erlangen) with a Flash EA 2000 Elemental Analyzer coupled to a Thermo Finnigan Delta V Plus mass spectrometer on samples with a spacing between ∼0.1 and ∼2.3 m. Measurement accuracy was ±0.06 ‰. Results are reported on the delta scale as δ¹⁸O, δ¹³C_carbonate and δ¹³C_org against the Vienna PeeDee Belemnite (VPDB) standard.

Initial stable isotope data from Cores LB1 and LB3 were presented in Lorenzen et al. (2013) and Moullade et al. (2015). Published data have a spacing of ∼20 cm in LB1 and ∼40 cm in LB3. New samples were analysed to increase the resolution to ∼5 cm during the early part of OAE1a between 18 and 41.5 m in Core LB1. Stable isotopes of bulk carbonate samples were analysed at the Leibniz Laboratory for Radiometric Dating and Stable Isotope Research at Kiel University with a Finnigan MAT 251 mass spectrometer. Analytical uncertainty for stable carbon isotopes was ±0.04 ‰ and ±0.07 ‰ for stable oxygen isotopes. Results are reported on the delta scale as δ¹⁸O and δ¹³C_carbonate against the Vienna PeeDee Belemnite (VPDB) standard.

2.5 Total organic carbon, carbonate and nitrogen content

Aliquots of samples for the determination of major and trace elements as well as phosphorus speciation in Core SN^∘4 were analysed for total organic carbon (TOC), carbonate and nitrogen content. Replicates of 8–10 mg sample material were pulverized, homogenized and sealed into tin capsules for the determination of total carbon (TC) and nitrogen (N) and into silver capsules followed by decarbonatization with 0.25 N hydrochloric acid (HCl) to measure TOC. Different weights of acetanilid (10.36 % N, 71.09 % C) and of the certified soil standard BSTD1 (0.216 % N, 3.5 % C) were sealed and measured for calibration and the monitoring of long-term stability. All capsules were measured with a Carlo Erba Protein Analyzer NA1500 at the GEOMAR Helmholtz Centre for Ocean Research Kiel. Carbonate content (TIC) was calculated as TIC = (TC-TOC) ×8.3333.

2.6 Time series analysis

The package astrochron (Meyers and Sageman, 2007; Meyers et al., 2012a; Meyers, 2014) for R (R Core Team, 2018) was used to perform frequency analysis on the NGR and XRF log(Terr∕Ca) datasets. The NGR datasets were linearly interpolated to a spatial resolution of 0.1 m for Core SN^∘4 and 0.025 m for the spliced dataset of Cores LB1–LB3. XRF-scanner-derived log(Terr∕Ca) from LB3 and LB1 was linearly interpolated to a spatial resolution of 0.01 m. Long-term variability of > 20 m wavelengths was removed from the NGR dataset by using a Gaussian kernel smoother with the function noKernel of astrochron.

Dominant frequencies in Cores LB1 and LB3 were extracted with the robust red noise multitaper method (MTM) spectral analysis of Mann and Lees (1996) using the function mtmML96 (Patterson et al., 2014) of astrochron with five 2π prolate tapers. Visualization and checks were performed with the eha function of astrochron with three 2π prolate tapers and a window of 6 m over the frequency range between 0 and 4 cycles per metre for the spliced record of LB1 and LB3. Identical intervals and frequency ranges were used to check results with the complementary XRF-scanner-derived log(Terr∕Ca) (Fig. S3.1).

An age model for Core SN^∘4 was derived by correlating the NGR record with that of the neighbouring Core S13 (Meyers et al., 2012a) following Kuhnt et al. (2017). A total of 14 tie points delimiting characteristic features of both cores were used (Sect. S4). The age model of Meyers et al. (2012a) was transferred into chronological ages by anchoring the Cenomanian–Turonian boundary at the top of cycle 3 to the chronological age of 93.9 Ma (Meyers et al., 2012b).

Periodicities for precession, obliquity and eccentricity are based on the orbital solution La04 (Laskar et al., 2004). Additional periodicities for eccentricity were extracted from the orbital solutions La10a-d (Laskar et al., 2011a) for the Cenomanian and Aptian intervals and La10a-d and La11 (Laskar et al. 2011a, b) for the Cenomanian interval using the astrochron function mtmML96 (Patterson et al., 2014). The analysis of orbital parameters for OAE1a was limited to the Aptian interval between 113 and 126 Ma and to the Cenomanian interval between 93 and 99 Ma for OAE2 (Sect. S5). Ratios were calculated for the main frequencies. Periodicities and ratios of the different orbital cycles are presented in Supplement Fig. S5.1 and Table S5.2 for the Aptian and in Table S5.3 for the Cenomanian interval.

2.7 Phosphorus speciation and analysis of major and trace elements

The distribution of major and trace elements as well as phosphorus speciation were determined on new samples from Core SN^∘4, mainly taken close to the published organic stable isotope samples (Beil et al., 2018). Samples cover the interval between 42.11 m (Section 18, Segment 3, 23–24 cm) and 305.77 m (Section 111, Segment 2, 29–31 cm). The core was sampled every ∼2.4 m with an increased resolution of ∼1.2 m over the global isotope excursions within OAE2 and the MCE. The surface of samples was cleaned with a metal-free porcelain knife and ground down with a rotary mill with agate balls to prevent metallic contamination. The powdered material was subdivided into aliquots for major and trace elements analysis and material for the determination of phosphorus speciation. A further aliquot was used for carbonate and organic carbon measurements for samples not close to the already published organic isotope samples. Additional samples from high-carbonate-content intervals were prepared and measured to determine the influence of low clay and organic content.

2.7.2 Phosphorus speciation

Samples for the measurement of phosphorus speciation were aliquots from samples used for the analysis of major and trace elements. A modified CONVEX extraction method (Oxmann et al., 2008; Sect. S6) was used to extract Ca-bound as well as Al∕Fe-oxyhydroxide-bound phosphorus. The sample material (250 mg) was subsequently treated with 3.75 mL KCl∕EtOH and decanted, and the supernatants were discarded three times. Al∕Fe-bound phosphorus was extracted from the centrifuged residue with the addition of 3.75 mL NaOH∕Na₂SO₄ (incubated for 1 h at 25 ^∘C), then 3.75 mL NaOH∕Na₂SO₄ (incubated for 2 h at 99 ^∘C) and finally 7 mL Na₂SO₄. Solutions were centrifuged and decanted after each of the three treatment steps and stored for measurement. The centrifuged residue was decarbonated for 8 h with 3.75 mL H₂SO₄. The Ca-bound phosphorus fraction was extracted by adding a further 3.75 mL H₂SO₄ (incubated for 2 h at 99 ^∘C) and twice Na₂SO₄ (3.75 and 7 mL). Each of the three steps was followed by centrifugation and decanting of the solution. Supernatants of each step were again mixed and stored for measurements. Details are provided in Fig. S6.1.

Figure 2Line scan photographs, stable isotope and XRF scanning elemental records spanning OAE2 in Tarfaya Basin Core SN^∘4. (a) Log(Terr∕Ca) from Kuhnt et al. (2017) and Beil et al. (2018). (b) Log(K∕Al) from Kuhnt et al. (2017) and Beil et al. (2018). (c) δ¹³O and (d) δ¹³C_org from Kuhnt et al. (2017) and Beil et al. (2018); a–c indicate prominent maxima following Voigt et al. (2007).

Download

For the measurement of Ca-bound phosphorus, an aliquot of 0.1 mL was diluted with 4.9 mL Milli-Q. Concentrations of Al∕Fe-bound phosphorus were measured on 4.5 mL of sample solution equilibrated with the addition of 0.5 mL of sol7 (Table S6.2) to a pH of ∼1. Concentrations of Al∕Fe- and Ca-bound phosphorus were calibrated with the PO₄ Merck Standard (1000 mg PO₄ L⁻¹) diluted to specific concentrations. The ammonium molybdate solution (0.1 mL; Table S6.3) was added 30 min before measurements. For samples close to the photometric saturation level a new aliquot was diluted to a lower concentration within the measurement range. A detailed list with the chemicals used and their respective concentrations is available in Table S6.2.

Figure 3Line scan photographs, stable isotope and XRF scanning elemental records spanning OAE1a from La Bédoule Cores LB3–LB1. (a) Log(Terr∕Ca). (b) Log(K∕Al). (c) δ¹³O and (d) δ¹³C_carbonate from Lorenzen et al. (2013) and Moullade et al. (2015) with new high-resolution isotope data from this study. The boundaries between individual phases (precursor, onset, peak, plateau, recovery) of OAE1a are indicated by dashed black lines. Carbon isotope stages C2–C8 of Menegatti et al. (1998) are indicated by red lines.

Download

The atomic ratio of organic carbon to P_reactive (defined as the sum of loosely absorbed, organic, authigenic and iron-bound phosphorus) was interpreted by Anderson et al. (2001) as a proxy for marine paleoproductivity. Due to sink switching (transfer to different phosphorus pools) between the different phosphorus species, especially during early diagenesis, the authors assumed that phosphorus in the different pools was originally derived from the remineralization of organic material and could therefore be considered equivalent to P_org. Here, P_react is calculated by summing up the two CONVEX-derived P pools: CaP (seen as equivalent to P_auth) and AlFeP. Both loosely absorbed and organic P were not measured. It is assumed that the concept of Andersen et al. (2001) is applicable, as both extracted P pools (CaP and AlFeP) contain the majority of P_total (mean =89 %; Fig. S7.1).

2.8 Phosphorus and organic carbon accumulation rates in Core SN^∘4

Bulk sediment mass accumulation rates (MARs) were calculated by multiplication of linear sedimentation rates (LSRs; derived from the age model of Core SN^∘4) with an average dry bulk density (DBD) of 2.1 g cm⁻³. The mean density for the interval encompassing OAE2 (600–2000 kyr) in neighbouring Core S13 (Meyers et al., 2012a) is 2.08 g cm⁻³ (SD 0.12 g cm⁻³). Phosphorus accumulation rates (PARs) were calculated by multiplying MAR with phosphorus concentrations from discrete measurements. Organic carbon accumulation rates (TOCARs) were calculated using the organic carbon concentrations from Kuhnt et al. (2017) and Beil et al. (2018).

3.1 Temporal evolution of bulk carbonate and organic carbon δ¹³C across OAE1a and OAE2

The high-resolution δ¹³C record of OAE1a in Core LB3–LB1 allowed for correlation with the Aptian isotope stages C2 to C8 proposed by Menegatti et al. (1998). We also identified prominent carbon isotope maxima a–c (following Voigt et al., 2007) in the SN^∘4 record of OAE2. We determined five phases (precursor, onset, peak, plateau and recovery) in the δ¹³C records of both OAE1a and OAE2 (Figs. 2, 3 and Sect. S8). The interval preceding both OAEs is characterized by a long-term mean without discernible trends and low variability. A precursor phase of two to four distinct δ¹³C minima (negative excursions) precedes the onset of the positive carbon isotope excursion (segment C3 of the Aptian carbon isotope curve of Menegatti et al., 1998). The onset phase encompasses the entire interval of increasing δ¹³C towards the first prominent peak (peak a of OAE2 following Voigt et al., 2007). This interval corresponds to segment C4 of the Aptian carbon isotope curve of Menegatti et al., 1998. The peak phase is defined as the interval starting with maximal values (first peak) to the end of the second peak (peak b of OAE2 following Voigt et al., 2007). A characteristic plateau phase (interval of relatively constant δ¹³C) follows the second peak of the carbon isotope excursions (C6 and lower part of C7 of Menegatti et al., 1998). The recovery phase encompasses the return of δ¹³C values to background levels following the OAE excursions (upper part of C7 and C8 of Menegatti et al., 1998).

The δ¹³C decreases during the precursor phase have comparable amplitudes, with 0.7 ‰ for OAE1a and 0.4 ‰ for OAE2 (Figs. 2 and 3). The most positive values of the OAE2 δ¹³C isotope excursion are reached at the first peak (peak a following Voigt et al., 2007) at 103.22 m with −24.2 ‰. By contrast, values continue to increase during OAE1a, with the highest values of 4.9 ‰ during the plateau stage at 40.11 m. The difference between the last minimum of the precursor phase and the first maximum (peak a for OAE2) of the peak phase is 2.4 ‰ for the carbonate δ¹³C record of OAE1a and 4.6 ‰ for the organic carbon δ¹³C record of OAE2. Considering the differing trend of OAE1a, with further increases during the plateau stage, the total amplitude of OAE1a is 3.6 ‰. Differences in mean values between the background levels before and after the OAE are 0.3 ‰ for OAE2 (from −28.5 ‰ (SD 0.3 ‰) before to −28.2 ‰ (SD 0.2 ‰) after OAE2) and 0.8 ‰ for OAE1a (from 1.9 ‰ (SD 0.1 ‰) before to 2.7 ‰ (SD 0.2 ‰) after OAE1a), indicating a general increase in background δ¹³C after the OAEs.

3.2 Temporal evolution of bulk carbonate δ¹⁸O across OAE1a and OAE2

The δ¹⁸O curves share common trends, except for the cyclic lithological changes in the upper part of the sedimentary record of OAE1a in LB3. Transient cooling events, identified by δ¹⁸O increases in LB3–LB1 and SN^∘4, occur during the early phases of both OAE1a and OAE2 (Figs. 2, 3 and Sect. S9). A first prominent cooling event prior to the onset of OAE2 in Core SN^∘4 (Kuhnt et al., 2017), which is not identified at other localities, was probably associated with local upwelling of cooler deepwater masses in the Tarfaya Basin. Cooling during OAE2 occurred in three main steps starting within the onset phase of the positive carbon isotope excursion. The most intense cooling, associated with the Plenus Cold Event, occurred during the peak phase of the excursion (in the trough between the δ¹³C peaks a and b). The Plenus Cold Event is globally recorded (e.g. Forster et al., 2007; Sinninghe Damsté et al., 2010; Jarvis et al., 2011; Jenkyns et al., 2017) and coincided with the invasion of boreal species in the European Chalk Sea (Gale and Christensen, 1996; Voigt et al., 2003), the extinction of the planktic foraminifer Rotalipora cushmani (e.g. Kuhnt et al., 2017) and reoxygenation of bottom water masses (e.g. Eicher and Worstell, 1970; Kuhnt et al., 2005; Friedrich et al., 2006). OAE1a shows a similar response of global temperatures to enhanced organic carbon burial (Kuhnt et al., 2011; Jenkyns, 2018): the main δ¹⁸O increase during the latter part of segment C4 in the δ¹³C curve of Meneghatti et al. (1998) also occurs during the onset phase. A further cooling event within segment C6 follows transient warming during the peak phase. Jenkyns (2018) recognized these transient coolings, also recorded in the northeastern Atlantic Ocean (Naafs and Pancost, 2016), Italy (Bottini et al., 2015), Turkey (Hu et al., 2012) and in the Pacific Ocean (Dumitrescu et al., 2006), as global events. Similarities in the δ¹⁸O records across both OAEs imply a similar response of the ocean–climate system to lowered atmospheric pCO₂ levels due to excess carbon drawdown associated with the burial of a vast amount of organic material on a global scale.

3.3 Regional differences in terrigenous input during OAE1a and OAE2

The different paleogeographic settings and depositional environments of the Tarfaya Basin (Core SN^∘4) and South Provence Basin (Cores LB1 and LB3) result in important differences in terrigenous sediment input. The XRF-derived log(Terr∕Ca) (Fig. 2) in the Tarfaya Basin exhibits low variability throughout, with the lowest values during OAE2. By contrast, log(Terr∕Ca) (Fig. 3) shows higher-amplitude variability during OAE1a. A major increase in log(Terr∕Ca) from −1.69 to −0.87 during the onset of OAE1a (68.09–67.99 m) indicates either a decrease in carbonate deposition or an increase in terrigenous input. Following this increase, the NGR and log(Terr∕Ca) records exhibit high-amplitude and high-frequency variability during C4 to C6 (26–23.29 m) and during the latter part of C7 and C8 (49.7–0 m) (Figs. 3, 4 and S3.1).

Table 2Durations and sedimentation rates for Aptian C stages (following Menegatti et al., 1998) in La Bédoule Core LB3–LB1 and comparison with durations from Malinverno et al. (2010) and Scott (2016).

Download Print Version | Download XLSX

3.4 Duration of OAE1a estimated from time series analysis of NGR and XRF scanner data

The spliced NGR record of LB1 and LB3 across OAE1a was subdivided into five intervals of relatively consistent orbital periodicities (N1 to N5) based on EHA (Fig. 4, Sect. S3). The recognition of cyclic patterns in the NGR record during OAE1a allowed for the correlation of major frequencies with orbital periodicities (Laskar et al., 2004, 2011a, b) and the calculation of mean sedimentation rates (Figs. 4 and S3.2, Tables S3.3 and S5.2). Comparison with durations proposed by Malinverno et al. (2010) and Scott (2016) for the different C stages of Menegatti et al. (1998) is shown in Table 2. A strong response to variations in obliquity and precession is evident in the latest part of stage C2, which precedes the onset of OAE1a (Figs. 4, S3.1 and S3.2). The NGR and log(Terr∕Ca) time series (Figs. 4 and S3.1) also exhibit precession-paced variations during the middle and later part of stages C7 and C8, which persist after OAE1a. This trend corresponds to pronounced lithological changes between carbonate-rich and clay-rich beds that are also apparent in the core images and the δ¹³C_carbonate and δ¹⁸O profiles (Fig. 3).

Figure 4Identification of orbital periodicities across OAE1a using evolutive harmonic analysis (EHA) of natural gamma ray (NGR) data in La Bédoule Cores LB3–LB1. (a) NGR; purple dashed lines mark boundaries from segment N1 to N4 (Sect. S3 and Fig. S3.2); mean sedimentation rates given in brackets. (b) EHA plot of NGR; orbital periodicities are indicated by white dashed lines; eccentricity e1 (404 kyr), e2 (124 kyr) and e3 (95 kyr), obliquity o1 (47 kyr) and o2 (37 kyr), precession p1 (22 kyr) and p2 (18 kyr) (Table S5.2). (c) δ¹⁸O and (d) δ¹³C_carbonate from Lorenzen et al. (2013) and Moullade et al. (2015) with new high-resolution isotope data from this study. Carbon isotope stages C2–C8 of Menegatti et al. (1998) are indicated by red dashed lines. Individual phases (precursor, onset, peak, plateau, recovery) of OAE1a are separated by black dashed lines. Note: the EHA colour scheme is limited to 50 % of the maximum amplitude to enhance visibility in the lower part of the core.

Download

Sedimentation rates (Table S3.3) increase from C3 and C4 (1.4 cm kyr⁻¹) to the recovery phase C8 (5.1 cm kyr⁻¹). The relatively stable ratio between terrigenous and carbonate content (log(Terr∕Ca) suggests a continuous increase in both carbonate and terrigenous input (Fig. 3).

3.5 Duration of OAE2 estimated from time series analysis of NGR data

EHA of NGR data reveals a marked response to obliquity (mainly o1, 48 kyr) and short eccentricity (e2 and e3; 100 kyr) during the positive δ¹³C excursion (Fig. 5). Following this interpretation, the δ¹³C excursion of OAE2 in Core SN^∘4 lasted for ∼790 kyr including the recovery phase (∼360 kyr without the recovery phase). The response to orbital forcing was more pronounced during the peak phase (∼90 kyr) and the recovery phase (∼440 kyr) than during the plateau phase (∼200 kyr). Sedimentation rates before the onset of OAE2 were 8.2 cm kyr⁻¹, dropping to 4.1 cm kyr⁻¹ during the positive δ¹³C excursion and recovering to ∼8 cm kyr⁻¹ in the upper part of the plateau phase (Fig. S4.1), which is comparable to levels prior to OAE2.

Table 3Durations of individual phases for OAE2 in Tarfaya Basin Core SN^∘4 and comparison with Li et al. (2017b) and Gangl et al. (2019).

Download Print Version | Download XLSX

Figure 5Identification of orbital periodicities over OAE2 from evolutive harmonic analysis (EHA) of natural gamma ray (NGR) data in Tarfaya Basin Core SN^∘4. (a) NGR. (b) EHA plot of NGR; orbital periodicities are indicated by white dashed lines: eccentricity e2 (127 kyr) and e3 (97 kyr), obliquity o1 (48 kyr) and o2 (38 kyr), precession p1 (22 kyr) and p2 (18 kyr) (Table S5.3). (c) δ¹⁸O and (d) δ¹³C_org from Kuhnt et al. (2017) and Beil et al. (2018); a–c indicate prominent maxima following Voigt et al. (2007).

Download

Four δ¹³C_org minima (precursor phase) occur over a period of ∼130 kyr prior to the positive δ¹³C excursion of OAE2 (Figs. 2, S4.1). The onset phase of OAE2 (∼70 kyr), defined as the interval between the centre of the last δ¹³C_org minimum and the first δ¹³C_org maximum of OAE2, is subdivided into three phases: the initial steep increase lasted ∼31 kyr, the intermediate plateau ∼14 kyr and the final rise to the first peak (peak a of Voigt et al., 2007) of the carbon isotope excursion ∼24 kyr. The time interval between the two maxima in δ¹³C_org (peak phase) was ∼90 kyr. The duration of the plateau phase, defined as the period between the second maximum and the end of the δ¹³C_org plateau, was ∼200 kyr. The recovery phase between the end of the plateau and the return to background values after OAE2 lasted ∼440 kyr. A comparison with durations proposed by Li et al. (2017b) and Gangl et al. (2019) is provided in Table 3.

3.6 Comparison of the durations of OAE1a and OAE2

The precursor phase of OAE1a (C3 of Menegatti et al., 1998), which consists of an extended interval of low δ¹³C_carbonate (< 1.75 ‰) and encompasses two prominent δ¹³C minima (< 1.5 ‰), lasted ∼434 kyr. The precursor phase of OAE2 is characterized by four nearly equally spaced δ¹³C_org minima (< 28.5 ‰) and lasted ∼126 kyr. The last two prominent minima immediately precede the onset of OAE2 and lasted for ∼75 kyr. The onset phase of the positive isotope excursions had a duration of ∼388 kyr for OAE1a (C4 of Menegatti et al., 1998) and ∼68 kyr for OAE2. The peak phase lasted ∼596 kyr in OAE1a (C5 and C6) and ∼86 kyr in OAE2, the plateau phase extended over ∼1343 kyr in OAE1a (C7) and ∼204 kyr in OAE2, and the final recovery phase had a duration of ∼377 kyr for OAE1a (C8) and ∼435 kyr for OAE2. The precursor, onset and peak phases were consistently ∼5 times longer in OAE1a than in OAE2, whereas the plateau phase was more than ∼7 times longer. The recovery phases had the same duration for both OAEs, taking into account uncertainties in defining change points at the end of the δ¹³C excursion.

3.7 Sedimentary phosphorus

Total phosphorus (P_total) concentrations in Core SN^∘4 (Fig. 6) are overall below 5 mg g⁻¹, except during two periods of enrichment peaking at 220.19 m with 16.61 mg g⁻¹ and at 104.4 m with 7.53 mg g⁻¹. Corresponding peaks in P_total ∕Al give a weight ratio of 1.07 at 220.19 m and 1.20 at 104.4 m. Concentrations and variability of P_total are low in the lowermost (305–267 m) and uppermost (67.51–42.11 m) parts of the studied interval, with means of 0.54 mg g⁻¹ (SD 0.19 mg g⁻¹) and 0.86 mg g⁻¹ (SD 0.3 mg g⁻¹), respectively. The origin of the phosphorus is difficult to assess, as the different species extracted are defined by the CONVEX method (Oxmann et al., 2008). Earlier studies found finely distributed fish debris and fecal pellets in Cenomanian sediments of the Tarfaya Basin (e.g. El Albani et al., 1999) partially reprecipitated as phosphate nodules during early diagenesis (e.g. Leine, 1986; Kuhnt et al., 1997). Phosphatic particles were not observed on the core surfaces and were not apparent in the XRF data from Core SN^∘4. They were only encountered as minor components in the residues of micropaleontological samples. The dissolution of phosphorus fixed in iron–aluminium crusts and the degradation of organic matter followed by reprecipitation as calcium-bound phosphorus further complicate the reconstruction of the initial source of phosphorus.

Figure 6Phosphorus concentration and speciation in Tarfaya Basin Core SN^∘4. (a) Concentration of total phosphorus (P_total). (b) Ratio of total phosphorus to aluminium (P∕Al) measured on bulk sediment. (c) Ratio of calcium- to aluminium- and iron-bound phosphorus (CaP∕AlFeP). (d) Concentration of aluminium- and iron-bound phosphorus (Al∕Fe-bound P). (e) Concentration of calcium-bound phosphorus (Ca-bound P). (f) Concentration of organic carbon (C_org) and (g) δ¹³C_org from Kuhnt et al. (2017) and Beil et al. (2018).

Download

3.8 Phosphorus speciation

Concentrations of reactive phosphorus (P_react; Fig. 7) in Core SN^∘4, calculated by summing up Al∕Fe- and Ca-bound phosphorus, are on average ∼89 % of the total phosphorus (P_total; Figs. 6 and S7.1) measured with ICP-OES. Organic-matter-bound phosphorus was not considered due to methodological limitations (Golterman, 2001) and to remineralization and reprecipitation (sink switching) into more stable Al∕Fe- and Ca-bound phosphorus species during early diagenesis. Differences between P_react and P_total are caused by non-extractable phosphorus bound in insoluble minerals or adhesively bound phosphorus extracted and discarded during the first step with KCl and EtOH (Fig. S6.1).

Figure 7Phosphorus speciation in Tarfaya Basin Core SN^∘4. (a) Concentration of reactive phosphorus (P_react) calculated as the sum of Al∕Fe- and Ca-bound P. (b) Concentration of organic carbon (C_org). (c) Ratio of aluminium- and iron-bound phosphorus to total phosphorus (AlFeP∕P_total). (d) Atomic ratio of organic carbon to reactive phosphorus (C_org∕P_react) with Redfield ratio 106:1 indicated as a dashed line. (e) Atomic ratio of total nitrogen to reactive phosphorus (N_total∕P_react) with Redfield ratio 16:1 indicated as a dashed line. (f) δ¹³C_org from Kuhnt et al. (2017) and Beil et al. (2018).

Download

Concentrations of phosphorus bound to aluminium or iron oxyhydroxides (AlFeP; Fig. 6) are low (median 0.014 mg g⁻¹). Increased concentrations (mean 0.027 mg g⁻¹) are determined for the lowermost interval between 305 and 272.58 m, followed by an interval with a lower average concentration (mean 0.014 mg g⁻¹) and lower variability (SD 0.006 mg g⁻¹) until 72.09 m. The uppermost interval is again characterized by increased variability (SD 0.008 mg g⁻¹). Calcium-bound phosphorus (CaP; Fig. 6) is the dominant species in the studied interval of Core SN^∘4. Concentrations are below 4 mg g⁻¹, except for two periods of enrichment peaking at 220.19 and 104.4 m with 15.76 and 6.62 mg g⁻¹, respectively. Both maxima coincide with the maximum enrichment of total phosphorus. The redox influence on Al∕Fe-bound phosphorus in Core SN^∘4 is addressed in the Sect. S10.

The atomic ratios of C_org and N_total against P_react show very similar trends except for the lower part between 305.77 and 258.35 m that is characterized in N_total∕P_react by increased values (median 5.9) and in C_org∕P_react by low values and low variability (median 77.6, SD 34.9), albeit in an interval with low resolution (Fig. 7). Afterwards both ratios show similar characteristics of low values (median of C_org∕P_react 89.1 and of N_total∕P_react 3.5) and low variability (SD of C_org∕P_react 67.1 and of N_total∕P_react 1.9) until 104.4 m, punctuated by short-lived maxima during the MCE (average of C_org∕P_react 170 and of N_total∕P_react 7.5). OAE2 is again characterized by markedly increased ratios (maxima of C_org∕P_react 1123.8 and of N_total∕P_react 29.6), and the uppermost part exhibits increased values (median of C_org∕P_react 242.6 and of N_total∕P_react 7.3). The ratio of CONVEX-extracted AlFeP to P_total (Fig. 7) shows increased values (mean 0.058) in the lower part between 305 and 272.58 m, followed by low ratios (mean 0.014) to the top of the sampled interval in Core SN^∘4.

3.9 Temporal changes in P-species concentrations and accumulation rates

Two maxima in P_total in the Cenomanian interval of Core SN^∘4 coincide with maxima in P_total ∕Al (Fig. 6). The more prominent of these phosphorus enrichments (peaking at 220.19 m) precedes the onset of the MCE by ∼3.04 m. The second maximum at 104.4 m occurs during the onset phase of OAE2 and coincides with a transient δ¹³C minimum, approximately halfway towards the first δ¹³C_org maximum at 103.22 m (Figs. 6 and S9.1).

The ratio of Ca-bound to Al∕Fe-bound phosphorus (Fig. 6) is characterized by three prominent peaks at 220.19, 143.51 and 104.40 m. The maxima at 220.19 and 104.40 m also exhibit high P_total concentrations and P_total ∕Al maxima. The first peak at 220.19 m precedes the first δ¹³C_org increase of the MCE, and the second peak at 104.4 m occurs within the onset phase of the δ¹³C_org excursion of OAE2. Both peaks of CaP∕AlFeP and P_total coincide with minima in the C_org∕P_total ratio (Fig. S7.1), suggesting an inorganic source for phosphorus or enhanced recycling of organic carbon and reprecipitation of the organic-bound phosphorus as Ca-bound phosphorus.

The increase in the atomic C_org∕P_total ratio (Fig. S7.1) during OAE2 post-dates the increase in the atomic C_org∕N_total ratio, interpreted by Beil et al. (2018) as enhanced cycling of nitrogen-rich organic matter within a dysoxic or anoxic water column. The increase in C_org∕P_total starts immediately above the prominent peaks in total and in Ca- and Al∕Fe-bound phosphorus, and it coincides with the first δ¹³C_org maximum of OAE2 (Figs. 6 and 7). The highest C_org∕P_total ratio is determined at 90.78 m within the plateau phase of OAE2. This peak coincides with minima in P_total and P_total∕Al, implying remobilization either synsedimentary due to preferential recycling of phosphorus in the water column and/or in the sediment. This decrease in P_total content is paralleled by a decrease in C_org (Fig. 6), suggesting enhanced remobilization of phosphorus as well as organic matter.

The atomic ratios of C_org∕P_react and N_total∕P_react (Fig. 7) are always lower than or close to the Redfield ratio ($\mathrm{C}:\mathrm{N}:\mathrm{P}$ $=\mathrm{106}:\mathrm{16}:\mathrm{1}$; Redfield, 1958, 1963), except during the MCE and OAE2, when C_org∕P_react is equal to or higher than 106:1 and N_total∕P_react is close to the predicted ratio of 16:1. Both ratios surpass predicted values at the first δ¹³C_org peak of OAE2, decrease slightly during the Plenus Cold Event and show a large increase during the plateau phase. The remaining interval is characterized by increased values above background level as in the lower part of Core SN^∘4 prior to the onset of OAE2.

Phosphorus accumulation rates (ARs) decline during the onset, peak and plateau phase of OAE2 (Fig. 8), with P_totalAR, P_reactAR and CaP AR declining by ∼75 % and AlFeP AR by ∼40 %. Phosphorus accumulation recovers after the end of the plateau phase and increases for P_total by ∼230 %, for P_react and CaP by ∼250 %, and for AlFeP by ∼210 %.

Figure 8Phosphorus accumulation rates across OAE2 in Tarfaya Basin Core SN^∘4 (age model based on Meyers et al. 2012a). (a) Accumulation rate of organic carbon (TOCAR). (b) Accumulation rate of aluminium- and iron-bound phosphorus (Al∕Fe-bound PAR). (c) Accumulation rate of calcium-bound phosphorus (Ca-bound PAR). (d) Accumulation rate of reactive phosphorus calculated as the sum of Al/Fe- and Ca-bound P (P_react AR). (e) Accumulation rate of total phosphorus (P_total AR). (f) δ¹³C_org. Stable isotope data and organic carbon concentrations from Kuhnt et al. (2017) and Beil et al. (2018); a–c indicate prominent maxima following Voigt et al. (2007).

Download

4.1 Influence of paleogeographic setting and weathering regime

Changes in the weathering regime of the source area during OAE1a, as shown by the XRF-scanner-derived log(K∕Al), influenced the sedimentary record of Cores LB1 and LB3 (Fig. 3). Increased log(K∕Al) and higher δ¹⁸O suggest predominantly physical weathering and/or intensified erosion, characteristic of drier conditions, and/or markedly seasonal rainfall prior to OAE1a (C2 and early C3) (Figs. 3 and 4). A decrease in log(K∕Al), synchronous with a shift to carbonate-depleted sediments indicated by log(Terr∕Ca), coincides with an increase in δ¹⁸O associated with the first transient cold event within the onset phase of OAE1a (C4). This major cooling event in the South Provence Basin was probably of global character (e.g. Jenkyns, 2018). The covariance of δ¹⁸O and log(K∕Al), in combination with published palynological and geochemical data (e.g. Masure et al., 1998; Hochuli et al., 1999; Keller et al., 2011; Föllmi, 2012; Cors et al., 2015), suggests that the La Bédoule area was located at the northern edge of the subtropical high-pressure desert belt and shifted into the Northern Hemisphere westerlies with increased rainfall during OAE1a. An equatorward contraction of the subtropical high-pressure belt also occurred in East Asia during the mid-Cretaceous warm period (Hasegawa et al., 2012). The paleo-location of La Bédoule (Cores LB1 and LB3) close to 30^∘ N (Masse et al., 2000) would predispose the area to latitudinal changes in the Hadley cell with drier conditions in times of an expanded high-pressure belt and increased rainfall resulting from an equatorward contracted Hadley cell. Strengthening fluctuations between drier and wetter climate conditions are recorded during the onset, main and early plateau phase. After a stabilization during the plateau phase, wet conditions prevailed after OAE1a.

The evolution of the weathering regime across OAE2 differed substantially in the Tarfaya Basin (Fig. 2), probably due to the lower paleolatitude and different paleoceanographic setting in a coastal upwelling area. The period prior to OAE2 was characterized by low sea surface temperatures and high mean log(K∕Al) exhibiting high variability, suggesting orbital forcing of a monsoonal hydrological cycle and weathering regime. A first prominent temperature decrease occurred during the late precursor and early onset phase of OAE2, similar to the global temperature decrease during OAE1a (Fig. 3; Jenkyns, 2018). A shift to lower mean log(K∕Al), indicative of stronger chemical weathering and/or reduced monsoonal seasonality in the source area, coincided with low sea surface temperatures during the precursor phase and suggests increased upwelling intensity and intensified monsoonal wind forcing. A shift back to higher log(K∕Al), indicating weaker chemical weathering in the source area, occurred at the end of the plateau phase of OAE2. The main interval of OAE2 (onset, peak and plateau phase) was characterized by intense chemical weathering and a weaker response of the hydrological cycle to orbital forcing, suggesting that the hinterland of the Tarfaya Basin was more or less permanently under the influence of tropical convective rainfall. The recovery phase and the immediate post-OAE2 period were once more characterized by weaker chemical weathering and a stronger response to orbital forcing, typical for a monsoonal regime at the northern edge of the seasonal swing of the Intertropical Convergence Zone (ITCZ). Persistently high temperatures (O'Brien et al., 2017) in the recovery and post-OAE2 phase, as indicated by low δ¹⁸O in SN^∘4, suggest that the equatorward shift of the monsoonal zone and southward expansion of the high-pressure desert belt occurred during a period of global warming. Changes in the weathering regime during the carbon isotope excursion may thus have been linked to major fluctuations in atmospheric CO₂ associated with increased carbon burial and enhanced response to local insolation forcing in the late onset and peak phase of OAE2 (Fig. 5).

4.2 Definitions and durations of OAE1a and OAE2

The estimated durations of OAE1a and OAE2 depend on the definitions of the events (based either on the stratigraphic extent of organic-rich sediments or δ¹³C excursions) and vary substantially from minimum estimates of 45 kyr for OAE1b to > 3 Myr for OAE1a (Table 1). Whereas there is a broader consensus in defining the duration of OAE2 as the entire interval between the onset of the positive δ¹³C excursion and the end of the recovery phase, the definition of OAE1a is commonly restricted to the interval of organic-rich sedimentation (Selli Level). This interval encompasses the onset and peak of the δ¹³C excursion but does not include the plateau and recovery phase (e.g. Malinverno et al., 2010). At La Bédoule, these differences result in durations of 1.1 Myr (onset and peak phases only) and 2.7 Myr (including plateau and recovery phases), which is almost 4 times the duration of OAE2. The differences in the definition of OAE1a and OAE2 stem from the different stratigraphic relationship between black shale accumulation and δ¹³C excursion in most of the “classic” localities in central and northern Italy, Britain, and North Africa. At these locations, the highest accumulation rates of organic carbon occurred within the peak and plateau phase of the carbon isotope excursion during OAE2, whereas during OAE1a the duration of black shale deposition was much shorter than the δ¹³C excursion.

4.3 Impact of orbital forcing on the evolution of OAE1a and OAE2

The orbital configuration favouring enhanced marine biological productivity in low latitudes may have been different for OAE1a and OAE2. Recent studies (e.g. Batenburg et al., 2016; Kuhnt et al., 2017) tuned the onset of OAE2 to a 405 kyr eccentricity maximum that succeeded an extended period of low seasonality, caused by a 2.4 Myr eccentricity minimum, which would be associated with stronger obliquity forcing during the precursor phase. The much shorter onset phase of OAE2 suggests a faster and stronger response of the ocean–climate system to the carbon cycle perturbation, either related to a shorter and more intense initial carbon dioxide release or to a different orbital configuration such as higher-amplitude eccentricity cycles. Additionally, higher rates of sea level rise would have promoted enhanced organic carbon burial in shallow shelf seas. These extensive dysoxic to anoxic shelf seas would have favoured the fast recycling of limiting nutrients, thereby enhancing primary production and further boosting organic carbon burial. Eccentricity control on both sea level and organic carbon burial is also suggested by the ∼90 and ∼600 kyr durations of the peak phases of OAE2 and OAE1a (maximum organic carbon deposition and the highest δ¹³C values) (Figs. 4 and 5).

By contrast, a strong obliquity imprint is detected during the plateau phases (Figs. 4 and 5). A period of high-amplitude obliquity forcing would have steepened the latitudinal temperature gradient and intensified atmospheric circulation (Batenburg et al., 2016; Kuhnt et al., 2017), thus helping to maintain elevated biological productivity in the tropical oceans. The obliquity-forced intensification of monsoonal systems may have resulted in periods of enhanced tropical weathering associated with enhanced nutrient supply to the ocean and wind-driven equatorial upwelling, which promoted carbon sequestration during the plateau and recovery phases. However, the recycling of essential nutrients from buried sediments or release from new sources (e.g. flood basalts of LIPs) would have been necessary to maintain this new equilibrium state over extended periods of time. The plateau stage lasted 7 times longer during OAE1a than during OAE2, which excludes orbital forcing as a primary control on its duration. This is also supported by differences in the response to orbital forcing at the end of the plateau phases: a precession and obliquity imprint during OAE1a contrast with a persistent obliquity signal during OAE2.

In addition to precessional variability, changes in orbital obliquity were a critical forcing factor of climate oscillations during the recovery phase of both OAEs. Obliquity determines the summer intertropical insolation gradient, which was recently suggested as an important driver of changes in tropical and subtropical hydrology and sedimentation patterns (Bosmans et al., 2015). An increased insolation gradient between the tropical summer and winter hemisphere during high obliquity leads to intensified atmospheric circulation within the Hadley cell, resulting in stronger cross-equatorial winds and intensified moisture transfer into the summer hemisphere. Continental climate proxy data and models suggest a contraction of the Hadley cell and a latitudinal shift of the subtropical high-pressure belt towards the Equator during mid-Cretaceous super-greenhouse conditions (Hasegawa et al., 2012; Hay and Flögel, 2012). This climate scenario would have placed the hinterland of the Tarfaya Basin in the dry, hot subtropical desert climate zone during OAE2, whereas the La Bédoule area was likely influenced by the humid zone of the Northern Hemisphere westerlies during OAE1a. An increase in cross-equatorial winds and intensified rainfall at obliquity maxima would thus have affected the Tarfaya Basin by increasing upper ocean mixing and upwelling and the La Bédoule area by increased rainfall due to intensified westerlies.

Periods of high-amplitude 41 kyr obliquity variations last ∼400–800 kyr and are separated by nodes of weak obliquity forcing with a duration of ∼200–400 kyr (Fig. S5.1, Laskar et al., 2004). During the recovery phase of OAEs, obliquity-forced low-latitude climate oscillations may have led to periods of enhanced equatorial upwelling, which would have promoted carbon sequestration and, over several hundred thousand years, depleted the ocean's nutrient pool. A period of low variability in orbital obliquity (obliquity node), commonly associated with global cooling episodes in Cenozoic climate records (e.g. Pälike et al., 2006), may have ultimately terminated OAE1a and 2.

4.4 Role of phosphorus recycling in maintaining high productivity during OAEs

Phosphorus is the primary limiting nutrient controlling marine biological productivity on longer (geological) timescales (e.g. Holland, 1978; Broecker and Peng, 1982; Smith, 1984; Codispoti, 1989), with the potential to control the occurrence of high-productivity events (e.g. Föllmi, 1996; Handoh and Lenton, 2003). In contrast to nitrate, which can be synthesized from atmospheric nitrogen primarily by cyanobacterial N₂ fixation under anoxic conditions (e.g. Rigby and Batts, 1986; Rau et al., 1987; Kuypers et al., 2004), the phosphorus supply to the ocean is restricted by riverine terrestrial input (Ruttenberg, 2003). This constitutes a limiting factor for increased marine primary productivity. Thus, alternative nutrient sources such as enhanced terrestrial input or the recycling of marine sediments were necessary to sustain high primary productivity over the extended durations of Cretaceous OAEs (e.g. Nedebragt et al., 2004).

C_org∕P_react and N_total∕P_react close to or above the Redfield ratio were measured within the MCE and OAE2 intervals, suggesting the preservation of the initial atomic ratio of primary production, possibly even enhanced by a change in the phytoplankton community (Geider and La Roche, 2002) or depletion in P_react. Unusually high phosphorus values in the intervals preceding the MCE and during the onset of OAE2 were likely caused by starved sedimentation, which may be due to the formation of hard grounds and condensed sections at sea level high stands or winnowing of fine-grained sediment by intensified bottom currents.

Regional changes in redox conditions towards anoxic or euxinic conditions in the lower water column would have inhibited the fixation of remobilized phosphorus in authigenic CaP, thus allowing the leakage of phosphorus into the water column and resulting in increased sedimentary C:P and N:P. Intensification of oxygen depletion during the plateau phase of OAE2 is associated with increasing C_org∕P_react and N_total∕P_react due to reduced diagenetic phosphorus precipitation, permitting phosphorus leakage from the sediment into the water column. However, phosphorus accumulation rates started to decline earlier during the precursor phase of OAE2. A trend toward more reducing conditions during the plateau and recovery phase of OAE2 was also suggested by Kolonic et al. (2005) based on high accumulation rates of redox-sensitive elements. A similar sequence of events was reconstructed by Stein et al. (2011) for OAE1a in deposits from the Gorgo a Cerbara section (Umbria–Marche basin) in central Italy (Sect. S11). Phosphorus accumulation and C_org∕P_total increased during the precursor phase (C3 of Menegatti et al., 1998) and decreased during the onset phase (C4 of Menegatti et al., 1998), suggesting enhanced phosphorus recycling from the sediments, as for OAE2 in Core SN^∘4 (Fig. 8 and Sect. S11). Stein et al. (2011) also reported a minor increase in C_org∕P_total at the base of the peak phase (C5+C6) and persistently high C_org∕P_total during the late peak phase (C5+C6), as for OAE2 in Core SN^∘4.

An earlier study by Poulton et al. (2015), focusing on Fe speciation proxies during the onset, peak and early plateau phase of OAE2 in nearby drill core S57, found cyclic variations between euxinic and ferruginous conditions. The stratigraphically extended interval (from the MCE to early Turonian) investigated at lower resolution by Scholz et al. (2019) is characterized by a high proportion of unpyritized reactive Fe in the total Fe pool. The results of this study are consistent with a proxy signature that is indicative of anoxic and non-sulfidic, so-called ferruginous, water column conditions throughout the studied interval (Poulton and Canfield, 2011). However, Scholz et al. (2019) argued that dissolved Fe (and hydrogen sulfide) concentrations in the water column of the Tarfaya system were unlikely higher than those observed in modern upwelling zones (e.g. Peru margin) on account of the low terrigenous sedimentation rates and tropical weathering on the adjacent continent. Previous studies proposed that phosphorus burial might be enhanced under ferruginous conditions, implying a negative feedback for the oceanic phosphorus pool and primary production (e.g. März et al., 2008). However, Scholz et al. (2019) did not observe a close relationship between Fe and P burial despite a ferruginous signature in the sediments, which supports the notion that dissolved Fe concentrations and rates of Fe oxide precipitation in the Tarfaya Basin were moderate and overall similar to modern upwelling systems (Wallmann et al., 2019).

The impact of oxygen depletion on the release of phosphorus into the water column was shown by recent studies on modern OMZs (e.g. Noffke et al., 2012; Schollar-Lomnitz et al., 2019) with implications for increased primary productivity through feedback mechanisms. Phosphorus recycling from the sediments may have sustained high primary productivity over extended periods of time, thus contributing to the long duration of OAE1a and OAE2. Phosphorus remobilization under anoxic conditions from vast areas of flooded shelf sediments facilitated high organic carbon burial rates during the onset, peak and plateau phases of the three anoxic events and may have acted as a positive feedback process, enhancing carbon burial and the removal of light carbon isotopes from the marine dissolved inorganic carbon reservoir and resulting in positive carbon isotope excursions.