We thank Anonymous Referee #1 for reviewing our manuscript. Here we reply to the comments given 14 Jun 2022, the revised manuscript will be changed accordingly.

First, we would like to point out, as noted by R#1, that For the first time, this high-quality and high time resolution dataset allows us to perform paleoceanographic studies on orbital time scale, probing the imprints of orbital forcing in the deep sea sediment of the Pacific Ocean during a warmer than today or greenhouse climate state. This is important because all sampling done on and observations obtained from the Shatsky Rise Late Cretaceous sediments are biased by aliasing. The new data we present for the first time clearly resolve Milankovitch cycles at the precession level. As such the key point of the manuscript is the precessional pacing observed as noted by R#1.

R#1 points out that looking at Figures S10, S12 and S13 the most prominent orbital cycle in the variability of planktonic d13C and d18O and the Ba is the 405 ka long eccentricity cycle, with strong 21 ka precession cycle, weak 41 ka obliquity cycle and relatively weak 100 ka short eccentricity cycle. And because the planktonic foraminiferal d13C is closely related to the biological pump and the carbon export production of the upper ocean R#1 asked why the manuscript deals with the precessional pacing of tropical ocean carbon export?

Based on the spectral analysis given in Figures S10, S12 and S13 the picture is more complex for each proxy data. Clearly, as given in the manuscript and used to establish the 405-kyr stable cyclostratigraphy, the bulk carbon stable isotope data are dominated by the long eccentricity cycle (highest peak in the MTM Fig. S10). A feature well known from the Cenozoic, very likely related to the long residence time of carbon in the ocean (See Pälike et al. 2006, Science). For the Late Cretaceous our record is the first that clearly shows the 405-kyr dominance in bulk carbon isotope data. The bulk oxygen isotope data (Fig. S12) do not show this strong 405-kyr influence but more power in the MTM spectra in the higher frequencies related to precession. One important point looking at the age MTM spectra is to be aware that the basic age model is a 405-kyr cyclostratigraphy which will enhance the long eccentricity related spectral peaks. Looking at the MTM spectra in the depth domain much less biased. However, there is a strong 405-kyr component in both carbon and oxygen bulk isotope data. We are cautious to interpret the oxygen bulk stable isotope data because it will be a mix of calcareous nannofossils from the surface ocean, planktonic foraminifera from the upper ocean layers and, to a lesser extent, the deep sea benthic foraminifera. It is not focus of this manuscript.

But the barium elemental data, which is focus of the manuscript, shown in the MTM spectrum in Figure S13 clearly is dominated by short cycles in the order of 20-30 cm and a really minor component of a 5m cycle, which is equivalent to the long eccentricity cycles if the short cycles are precession (~20 precession cycles in one 405-kyr cycle; ~20 25cm cycles in a 5 m cycle). Even with the 405-kyr age model the MTM spectrum clearly shows the domination of precession. Looking at Figure 3 of the main text (this Figure will be full page in the revised version to better illustrated the fine cyclicity) and the supplementary Figure S4 as well as the correlation Figure S8 there is one dominating rhythm related to precession with very little expression of modulations by eccentricity. Prominent examples of eccentricity modulated precession cycles in geochemical data is the late Paleocene and early Eocene (see Lourens et al. 2005; Westerhold et al. 2007 and 2008) are published from Walvis Ridge in the South Atlantic. Compared to those records the Maastrichtian Shatsky Rise XRF records show a minor eccentricity component. Because of the modulation of the precession by eccentricity some related variation in amplitude of the data can be expected and is seen, but not to an extent that would lead to the interpretation of major importance on the data.

In the manuscript we focus on the cyclicity of XRF Barium data which we show are very likely related to changes in surface ocean productivity. And because these data are dominated by precession our manuscript focuses on the precessional pacing of tropical ocean carbon export. We think this is justified by the outstanding quality of the data. No changes will be made to the manuscript following this comment.

The traditional concept of “biological pump” was cited to interpret the ocean carbon export production in this study. In Figure S5, the authors compared the high time resolution bulk d13C and d18O with the low time resolution planktonic foraminiferal d13C and d18O. They concluded that the two kinds of d13C records resemble to each other. As far as I am concerned, only the bulk carbonate d13C shows similarity to the Rugoglobigerina rugosa (subsurface species, deeper) d13C (Figure S5, top, green) but big differences with the Pseudoguembelina costulata and P. kempensis (surface species, shallower) d13C (Figure S5, top, orange). In the greenhouse world of the Cretaceous, a time interval without global ice volume effect on the calcareous shells of planktonic foraminifers, the temperature and salinity are two major factors controlling the planktonic foraminiferal d18O. As seen in Figure S5, the d18O of the subsurface species (green) are obviously heavier than that of the surface species (orange), indicating cooler or saltier water mass. However, the d13C of the subsurface species (green) is also heavier than that of the surface species (orange), which is opposite to the vertical distribution of the water mass d13C caused by the traditional biological pump that decreases with water depth. The inconsistency probably indicates that the ocean carbon export production of the tropical Pacific Ocean during the greenhouse world of the late Cretaceous is different with that we have known today. The authors need to add a new paragraph to discuss this inconsistency before they could use the XRF core scanning Ba record as a proxy of the ocean carbon export production.

Explaining this apparent inconsistency that the reviewer outlines in the foraminiferal isotope data does not preclude using the XRF scanning Ba record as a proxy of ocean carbon export production since the Ba proxy is independent of the foraminiferal and bulk C isotope data being that it relies on the modern observations that the mineral barite accumulates under regions of high organic carbon export can be used to independently reconstruct changes in organic carbon export following previous work in the Cenozoic (see Dymond et al., 1992; Paytan et al., 1996; Eagle et al., 2003; Griffith et al., 2021). The XRF scanning Ba record is confirmed in our study as a proxy for the mineral barite by carefully separating barite from discrete sediment samples (see Figs. 2, S9, S11). This is necessary to be sure that the XRF scanning record is actually recording marine pelagic barite and not a diagenetic phase or some other Ba-phase.

The relevance of the foraminiferal isotope records is in better understanding the dynamics in the water column – as the reviewer points out – and the primary influence on the bulk carbonate d13C record presented in high resolution for the first time in this study. Whether or not Rugoglobigerina rugosa represents surface or subsurface waters is less clear. Frank et al. (2005) suggested that R. rugosa occupied the mixed layer depth and contained photosymbionts (Abromovich et al., 2003)  and was the most enriched in 13C followed by two other planktic foraminifera measured (G. stuartiformis and P. multicamerata) – similar to our presentation of R. rugosa with P. costulata and P. kempensis all of which have been suggested as mixed layer dwellers (see Abramovich et al., 2003). This suggests that P. costulata and P. kempensis could have been somewhat deeper dwelling than R. rugosa at this time with a surface d13C gradient consistent with expectations of the heaviest d13C at the shallowest depths. Jung et al. (2013) suggested that their planktic R. rugosa d18O record (plotted in Fig. S5) might have been altered toward cooler temperatures during early diagenetic recrystallization taking place on the seafloor – or as Reviewer 1 suggested R. rugosa could be recording a saltier water mass.

We thought it was important to see if the d13C variations we measured in the bulk carbonate record at high resolution could be related to surface production – and found that they were generally similar to the planktic R. rugosa d13C record suggesting that they could be recording changes in local surface productivity. So we included the following:

Original manuscript (Lines 188-190): “The bulk carbonate δ13C values from Shatsky Rise closely resemble the lower resolution surface planktic foraminiferal δ13C records from the same sites reported by Jung et al. (2013) and Dameron et al. (2017). Hence, the bulk carbonate isotope record in the region primarily reflects surface conditions, which can be influenced by changes in local productivity (Fig. S5).”

To further clarify in the main manuscript what is seen in fig. S5, we will revise this paragraph to read in the revised manuscript: “The bulk carbonate δ13C values from Shatsky Rise closely resemble the lower resolution surface planktic foraminiferal δ13C record of Rugoglobigerina rugosa from the same sites reported by Jung et al. (2013). Hence, the bulk carbonate isotope record in the region likely reflects surface conditions, which can be influenced by changes in local productivity (Fig. S5).”

The figure caption for the supplement will outline what is mentioned above:

“Figure S5. (from top to bottom) Shatsky Rise composite bulk δ13C (black line) with 50 points moving average (red line) shown with planktic foraminifera δ13C from Rugoglobigerina rugosa (Jung et al., 2013; green circles) and planktic foraminifera δ13C from Pseudoguembelina costulata and P. kempensis (Clark et al., 2012; Dameron et al., 2017; orange circles). All three were thought to have occupied the mixed layer depth (Abromovich et al., 2003), however given that P. costulata and P. kempensis have lower δ13C values, they likely were somewhat deeper dwelling than R. rugosa at this time if the surface δ13C gradient was consistent with expectations of the heaviest d13C at the shallowest depths. Shatsky Rise composite bulk δ18O (black line) with planktic foraminifera δ18O. XRF Ba (total counts) with 100 points moving average (red line). At the bottom, XRF Sr/total area is plotted with planktic:benthic foraminifera ratio (P:B ratio; Dameron et al., 2017) in green. Note age is decreasing from the left to right.”

New reference will be added to Supplement: Abramovich, S., G. Keller, D. Stuben, and Z. Berner (2003), Characterization of late Campanian and Maastrichtian planktonic foraminiferal depth habitats and vital activities based on stable isotopes, Palaeogeogr. Palaeoclimatol. Palaeoecol., 202, 1–29.

The last part of this manuscript “3.4” focuses on the discussion of direct response to precession during greenhouse world. An important content of this part is the discussion on whether the tropical Pacific was more like a permanent El Niño like state or robust ENSO- like variability existed in past greenhouse conditions. ENSO is El Niño-southern oscillation. The authors need to explain the difference in “El Niño like state” and “robust ENSO variability” in a greenhouse condition. Today, we usually use the changes of the gradients in both the SST (Sea Surface Temperature) and the Thermocline Depth of the east and west equatorial Pacific to depict the ENSO variability that is a typical climate phenomenon in the equatorial Pacific Ocean. The sites 1209 and 1210 were in the middle of the tropical open ocean in the late Cretaceous (Figure 1). Can we depict the ENSO variability if it really existed in the late Cretaceous without reconstructions of the gradients of the SST or thermocline depth in the east and west equatorial ocean? If not, the vague discussion based on non-proxy-derived discussion would lead to misunderstanding.

Yes, we agree with the reviewer that in order to reconstruct ENSO variability at this time, more than one record is needed. However, we argue that if there was a permanent El Nino like state at this time in the Pacific (with deeper thermocline and more stratified tropical Pacific) as previously hypothesized for past warm climates, we would not see the large variations in export production at our tropical Pacific site on relatively short timescales, i.e., precessional changes of more than 2x in export production and more than 4x in benthic foraminiferal accumulation rates.  

Original manuscript with proposed revisions “…Previous work suggested a link between past warm climates and a “permanent El Nino state” (Wara et al., 2005; Fedorov et al., 2006) with a deeper thermocline and more stratified tropical Pacific preventing upwelling of nutrients that result in reduced primary production and low carbon export. However, the large range of variations in our carbon export record in the tropical Pacific (XRF-Ba content and benthic foraminiferal accumulation rates) during the Maastrichtian greenhouse argues against this hypothesis of a continual El Nino-like state. We suggest instead that this new record adds to the growing body of work that suggests robust El Nino-Southern Oscillation or ENSO-like variability existed in past greenhouse conditions (e.g., Davies et al., 2012).”

Add new reference to main manuscript: Wara, M.W., Ravelo, A.C., and Delaney, M.L., 2005, Permanent El Niño–like conditions during the Pliocene warm period: Science, v. 309, p. 758–761, doi:10.1126/science.1112596.

We thank Referee #2 Mingsong Li for reviewing our manuscript. Here we reply to the comments given 26 Jun 2022, the revised manuscript will be changed accordingly.

Two major and some minor issues were raised by R#2. First we would like to comment on the two major issues:

1) Missed information about the tuning in the main paper

R#2 suggests to insert a chapter about the astrochronology / cyclostratigraphy in the main text as this is a central part of the study. In a revised manuscript we will add a chapter dealing with the age model development. We will move supplementary figures S2 (Basic 405 kyr age mode), S10 (bulk carbon isotope MTM power spectrum and wavelet analysis in depth and age) and S12 (Barium MTM power spectrum and wavelet analysis in depth and age) to the main text.

R#2 points out that the astronomical tuning strategy is generally straightforward in this case study and that the resulting mean sedimentation rate of about 1.25 cm/kyr from the 405-kyr cyclostratigraphy is supported by other geologic evidence.

We have tested our 405-kyr cyclostratigraphy by plotting the Shatsky Rise bulk δ13C record against the bulk δ13C record of the Zumaia succession from Batenburg et al. (2012) which is the only record with comparable resolution. As given in the main text we updated the tie points of Batenburg et al. (2012) to the Laskar cosine function (Table S2) to make the records consistent with respect to the target curve. This way we can check if the number of 405-kyr cycles identified in the Shatsky Rise record is consistent with the complete record from Zumaia.

R#2 comments that the method we used can be subjective and it is unclear whether all assigned 405 kyr cycles are fine or not. One solution is to use complex statistical tools like TimeOpt or COCO, where the latter was developed by R#2. And indeed as given by the test by R#2 using COCO the results are basically the same as ours. We thank R#2 for point out to use these methods. Moreover, R#2 suggests a potential critical interval identified by his method at ~285m rmcd. This corresponds to between 405-kyr cycles labeled Maa5 and Maa6 in Figure S2 and S3. If there is one more 405-kyr cycle in this interval the carbon isotope minimum at 68.3 Ma would not match the Zumaia record anymore.

In the revised version of the manuscript we will add the COCO approach for testing our age model. It should be pointed out here that the drill cores used are affected by coring disturbance to some extent, particularly when chert layers were encountered. For the composite record from Shatsky Rise we tried to avoid those intervals whenever possible. However, there will be some disturbance and therefore misinterpretation of changes in cycle thickness can occur using an automated analysis routine that assumes constant sedimentation rates. Unpublished precession cycle by precession cycle correlation to Zumaia also reconfirms the basic 405-kyr age model presented in the initial manuscript.

1. Unclear relationship between the thermocline and precession

This paper argues that precession paced ocean export production via influencing the thermocline depth. In the text, shoaling of the thermocline is thought to have occurred in response to changes in precession. This is a critical component of the story; however, the mechanism is not presented. In Figure 5, the authors note that “a shallower thermocline at high precession” represents a huge jump. Please fulfill this gap in the main text.

Moreover, the paper mentioned that “ENSO-like variability” existed in past greenhouse conditions. It would be better if this variability pattern can be introduced so readers don’t have to read other papers (e.g., Davies et al., 2012) to understand this term.

We will add additional introduction to this variability in the revised manuscript and more clearly lay out the proposed relationship between precession and this ENSO-like variability (proposed revisions italicized)  “…Previous work suggested a link between past warm climates and a “permanent El Niño state” (Wara et al., 2005; Fedorov et al., 2006) with weak trade winds along the equator, a deeper thermocline and more stratified tropical Pacific preventing upwelling of nutrients that result in reduced primary production and low carbon export. Today this occurs as an irregular periodic variation along with an opposite state in the tropical Pacific with strong trade winds and more intense upwelling in the eastern tropical Pacific driving increases in productivity on interannual-to-decadal timescales. Changes in the mean state over longer timescales can be recorded in deep sea sediment archives (e.g., Pena et al., 2008; Zhang et al., 2021). The large range of variations in our carbon export record in the tropical Pacific (XRF-Ba content and benthic foraminiferal accumulation rates) during the Maastrichtian greenhouse argues against the hypothesis of a continual El Niño-like state. The new record exhibits variations which we suggest reflect changes in the mean climate state. We suggest instead that this new record adds to the growing body of work that suggests robust El Niño-Southern Oscillation-like variability existed in past greenhouse conditions (e.g., Davies et al., 2012) and may be sensitive to orbital forcing, especially the effect of orbital precession in the tropics (e.g., Clement et al., 1999; 2000; Lu et al., 2019). The resultant tropical mean climate state and climate variability forced by minima and maxima of precession changing the strength of coupled ocean-atmosphere feedbacks in the tropical Pacific is less clear during the Maastrichtian greenhouse (Fig. 4), and requires focused modeling efforts and new proxy records which test these hypothesized relationships.”

Upon further examination and reflection of the literature, we want to shift the focus of the discussion here (see above). Figure 4 will be revised so that a precession is more generally related to changes in the mean state of the system, shifting from a mean state with a shallower thermocline resulting in higher productivity to a deeper thermocline with lower productivity in the region. Without additional work (outside the scope of this study) it is not possible to know if it is minima or maxima in precession that should result in one mean state or the other.

New references to be added to the main manuscript:

Clement, A. C., Seager, R., and Cane, M. A.: Orbital controls on the El Niño/Southern Oscillation and the tropical climate, Paleoceanography, 14, 441– 456, 1999. doi:10.1029/1999PA900013.

Clement, A. C., Seager, R., and Cane, M. A.: Suppression of El Niño during the mid-Holocene by changes in the Earth's orbit, Paleoceanography, 15, 731– 737, 2000. doi:10.1029/1999PA000466.

Liu, Z., Lu, Z., Wen, X., Otto-Bliesner, B., Timmermann, A., and Cobb, K.: Evolution and forcing mechanisms of El Nino over the past 21,000 years, Nature, 515(7528), 550– 553, (2014). https://doi.org/10.1038/nature13963

Lu, Z., Liu, Z., Chen, G., and Guan, J.: Prominent precession band variance in ENSO intensity over the last 300,000 years, Geophysical Research Letters, 46, 9786-9795, 2019. doi:10.1029/2019GL083410.

Pena, L. D., Cacho, I., Ferretti, P., and Hall, M. A.: El Niño-Southern Oscillation-like variability during glacial terminations and interlatitudinal teleconnections, Paleoceanography, 23, PA3101, 2008. doi:10.1029/2008PA001620

Wara, M. W., Ravelo, A. C., and Delaney, M. L.: Permanent El Niño–like conditions during the Pliocene warm period, Science, 309, 758–761, 2005. doi:10.1126/science.1112596.

Zhang, S., Yu, Z., Gong, X., Wang, Y., Chang, F., Lohmman, G., Qi, Y., and Li, T.: Precession cycles of the El Niño/Southern oscillation-like system controlled by Pacific upper-ocean stratification, Communications Earth & Environment, 2, 239, 2021. doi:10.1038/s43247-021-00305-05.

Minor issues

1. Line 28: add “405 kyr,” before 100 kyr.

Will be corrected in the revised manuscript

1. Line 111: sampling time of 12 s. Do you have evidence that 12 s is sufficient for the data quality? I have tested a range of sampling times ranging from 10 s to 120 s, it looks like XRF can usually produce stable results after ca. 20 s.

The Shatsky Rise cores consist of sediments with very high carbonate content, typically >95%. XRF scanning of very high carbonate sediments is difficult due to artifacts produced by the extremely high calcium peak and interaction with other elements. Thus, after test scans, T. Westerhold decided to scan only at a 50 kV setting because only the barium signal gave a stable and good result using 12 seconds count time. Agreed that count time should be longer to get statistically better results, however there is a time and cost limitation running cores at the repositories. Thus 12 seconds only scanning 50 kV was deemed useful to establish a stratigraphy and look at productivity related changes. Other elements like iron, which normally are pretty robust in acquiring, are difficult to measure on the high carbonate cores in relatively short time. One has to be aware that >95% carbonate content at a 10 kV run (usually the one to measure Al to Fe) will result is a high deadtime prolonging the scan time enormously. From the experience of T. Westerhold with many XRF scanned sections on a large range of sediment types and cores, stable results can be obtained already at 5-7 seconds count time having the right sediment composition and sediment core quality (not too high water content etc.). Generally, it cannot be stated as done by R#2 that XRF can usually produce stable results after ca. 20 s only.

1. Figure 1: Please explain the meaning of shade zones in the map (brown, light olive, blue, and purple).

Will be included in the revised manuscript

1. Figure 3: obliquity period was not 41 kyr. Alternatively, the La2004 solution predicts a 38 kyr cyclicity for this time interval.

Will be corrected in the revised manuscript

1. Figure 4: the middle rows are too busy. Why there is no BFAR data at the interval of 66.74- 68.71 Ma?

BFAR was only done for two cycles (peaks in Ba) to confirm whether or not the organic matter export indicated by the Ba peaks was reaching/impacting the seafloor community. No data was collected for their earlier peak.

1. Line 468: Remove “Strong” in the caption.

Will be corrected in the revised manuscript