Paleogeographic controls on the evolution of Late Cretaceous ocean circulation

Understanding of the role of ocean circulation on climate during the Late Cretaceous is contingent on the ability to reconstruct its modes and evolution. Geochemical proxies used to infer modes of past circulation provide conflicting interpretations for the reorganization of the ocean circulation through the Late Cretaceous. Here, we present climate model simulations of the Cenomanian (100.5–93.9 Ma) and Maastrichtian (72.1–66.1 Ma) stages of the Cretaceous with the CCSM4 earth system model. We focus on intermediate (500–1500 m) and deep (> 1500 m) ocean circulation and show that while there is continuous deep-water production in the southwestern Pacific, major circulation changes occur between the Cenomanian and Maastrichtian. Opening of the Atlantic and Southern Ocean, in particular, drives a transition from a mostly zonal circulation to enhanced meridional exchange. Using additional experiments to test the effect of deepening of major ocean gateways in the Maastrichtian, we demonstrate that the geometry of these gateways likely had a considerable impact on ocean circulation. We further compare simulated circulation results with compilations of ε_Nd records and show that simulated changes in Late Cretaceous ocean circulation are reasonably consistent with proxy-based inferences. In our simulations, consistency with the geologic history of major ocean gateways and absence of shift in areas of deep-water formation suggest that Late Cretaceous trends in ε_Nd values in the Atlantic and southern Indian oceans were caused by the subsidence of volcanic provinces and opening of the Atlantic and Southern oceans rather than changes in deep-water formation areas and/or reversal of deep-water fluxes. However, the complexity in interpreting Late Cretaceous ε_Nd values underscores the need for new records as well as specific ε_Nd modeling to better discriminate between the various plausible theories of ocean circulation change during this period.

1 Introduction

Over the last several decades, a wealth of proxy data established that the Cretaceous period was characterized by a greenhouse climate, with warmer-than-modern temperatures and an absence of perennial polar ice sheets (e.g., Barron, 1983; Jenkyns et al., 2004; O'Brien et al., 2017). This characterization draws on paleontological and paleobotanical data, including the findings of fossils of ectothermic species (e.g., Tarduno et al., 1998) and woody vegetation (e.g., Bowman et al., 2014) at polar latitudes, as well as geochemical studies indicating warm sea surface and deep ocean temperatures at all latitudes (e.g., Wilson and Norris, 2001; Pucéat et al., 2003; Friedrich et al., 2012; MacLeod et al., 2013; O'Brien et al., 2017; Huber et al., 2018). Successive refinements of the data indicating Cretaceous warmth also reveal a greater variability within Cretaceous climates than often portrayed and includes carbon cycle perturbations referred to as ocean anoxic events (OAEs; e.g., Schlanger and Jenkyns, 1976; Jenkyns, 2010) and intervals of cooler climatic conditions indicated by evidence for polar sea ice (Davies et al., 2009; Bowman et al., 2013) and possibly short-lived polar ice sheets (Price, 1999; Ladant and Donnadieu, 2016). Global paleotemperature compilations confirm this variability and provide evidence of global warming through the Early Cretaceous to early Late Cretaceous (Cenomanian–Turonian) interval of maximum temperatures followed by cooling through the end of the Cretaceous (Cramer et al., 2011; O'Brien et al., 2017; Huber et al., 2018).

Early attempts at modeling past climates with atmosphere-only global climate models suggested that Cretaceous warmth was the result of paleogeographic changes and higher atmospheric CO₂ concentrations (Barron and Washington, 1982, 1984, 1985). The role of paleogeographic changes on global temperature evolution across the Cretaceous has been debated for a long time (Poulsen et al., 2001; Donnadieu et al., 2006; Fluteau et al., 2007). Recent experiments with models of higher complexity and higher-resolution support only a weak impact of Cretaceous paleogeographic changes on global temperature evolution (Lunt et al., 2016; Tabor et al., 2016). In contrast, model simulations (Poulsen and Zhou, 2013; Tabor et al., 2016) and paleo-CO₂ reconstructions (Fletcher et al., 2008; Wang et al., 2014) suggest that atmospheric CO₂ concentrations provided a first-order control on Late Cretaceous temperatures. Indeed, compilations of paleo-CO₂ concentrations across the Cretaceous suggest that CO₂ and temperatures broadly increased to peak levels during the Cenomanian and Turonian thermal maximum, before decreasing throughout the Late Cretaceous (Wang et al., 2014). The comparison between model simulations and proxy reconstructions of sea surface temperatures (SST) provides further support for a Late Cretaceous cooling trend driven by decreasing CO₂ levels (Tabor et al., 2016).

Late Cretaceous cooling is expressed heterogeneously at a regional scale and reveals interbasinal variations in both the surface and deep ocean (Friedrich et al., 2012; O'Brien et al., 2017; Huber et al., 2018). For instance, records from the North Atlantic and Indian Ocean show cooling from the Turonian to the mid-Campanian and stabilization or warming afterward, whereas records from the Pacific Ocean and from the Atlantic and Indian sectors of the Southern Ocean show gradual cooling from the Turonian to the Maastrichtian (e.g., MacLeod et al., 2005; Huber et al., 2018). These distinct regional trends suggest that the pathways followed by water masses and connections between ocean basins changed during the Late Cretaceous as a result of the evolving paleogeography.

This conjecture is corroborated by studies of the temporal trends and spatial variations in the neodymium (Nd) composition of seawater (i.e., the ratio of ¹⁴³Nd ∕ ¹⁴⁴Nd, expressed as ε_Nd). Seawater ε_Nd values are mainly controlled by export of dissolved Nd through continental weathering and fluvial runoff to the ocean (e.g., Frank, 2002; Goldstein and Hemming, 2003; Tachikawa et al., 2017), but mass-balance calculations have shown that additional sources, such as exchange with continental margins (or boundary exchange, Lacan and Jeandel, 2005), are required to close the Nd budget (Tachikawa et al., 2003; Lacan and Jeandel, 2005; Arsouze et al., 2007; Rempfer et al., 2011). Because the residence time of Nd in the ocean is shorter than the oceanic mixing time, and initial Nd isotopic ratios are not totally overprinted by particle and/or boundary exchange during circulation (e.g., Frank, 2002; Tachikawa et al., 2003; Rempfer et al., 2011), the ε_Nd composition of water masses acts as a quasi-conservative tracer reflecting the geographical provenance and oceanic pathway of water masses (Piepgras and Wasserburg, 1982; Frank, 2002; Goldstein and Hemming, 2003; Moiroud et al., 2016; van de Flierdt et al., 2016) and, as such, is used as a proxy for past ocean circulation.

Records of Nd isotopes illustrate, in particular, a long-term shift toward more unradiogenic (lower) values in the Atlantic basin between the Turonian and the Campanian (e.g., Robinson et al., 2010; MacLeod et al., 2011; Martin et al., 2012; Robinson and Vance, 2012; Moiroud et al., 2016; Batenburg et al., 2018). However, there is no consensus on the specific modes and evolution of ocean circulation across the Late Cretaceous as interpretation is complicated by the lack of Late Cretaceous ε_Nd records in key places and times, by the possibility of modification of ε_Nd values along flow paths, and by uncertainties in the paleodepth of sites where ε_Nd values were documented. Illustrating this lack of consensus, deep-water formation during the Late Cretaceous has been hypothesized to occur (alternatively or concurrently) in most high-latitude basins, including the South Atlantic and Indian Ocean (e.g., Robinson et al., 2010; Murphy and Thomas, 2012; Robinson and Vance, 2012), North Atlantic (e.g., MacLeod et al., 2011; Martin et al., 2012), North Pacific (e.g., Hague et al., 2012; Thomas et al., 2014; Dameron et al., 2017), and South Pacific (e.g., Thomas et al., 2014; Dameron et al., 2017; Haynes et al., 2020), as well as possibly in the low latitudes (e.g., Friedrich et al., 2008; MacLeod et al., 2008, 2011).

Numerical models have been instrumental in providing a framework for interpreting the paleoceanographic data and in shedding light on new hypotheses, yet the location of possible sources of deep water differs among simulations almost as much as it does among conclusions of proxy studies (e.g., Brady et al., 1998; Poulsen et al., 2001; Otto-Bliesner et al., 2002; Zhou et al., 2008; Monteiro et al., 2012; Donnadieu et al., 2016; Ladant and Donnadieu, 2016; Lunt et al., 2016). Numerous factors may explain this spread, in particular differences in model complexity and resolution and differences in the paleogeography employed, which may vary across model studies (Donnadieu et al., 2016; Lunt et al., 2016; Tabor et al., 2016). Even within identical Cretaceous time slices, there can be significant differences in paleogeographic reconstructions resulting in additional uncertainty regarding the areas of deep-water formation as well as the configuration of oceanic gateways and, thereby, the modes of ocean circulation (e.g., Donnadieu et al., 2016; Lunt et al., 2016; Farnsworth et al., 2019). The considerable impact of breaching a continental barrier or closing an oceanic seaway has long been demonstrated in idealized and paleoclimatic model studies (e.g., Toggweiler and Samuels, 1995; Poulsen et al., 2003; Sijp and England, 2004; Sepulchre et al., 2014; Donnadieu et al., 2016; Elsworth et al., 2017; Ladant et al., 2018; Tabor et al., 2019).

Interbasinal differences in temperature evolution and shifts in the global ocean circulation, therefore, point toward a critical role of paleogeographic reorganizations, such as the geometry of oceanic basins or the opening, closure, and depth changes in oceanic gateways, regardless of the evolution of atmospheric CO₂ during the Late Cretaceous. To our knowledge, only one coupled ocean–atmosphere model study focused on the evolution of global ocean circulation during the Late Cretaceous (Donnadieu et al., 2016). Using Cenomanian–Turonian and Maastrichtian simulations, Donnadieu et al. (2016) demonstrated a shift toward a more vigorous ocean circulation in the Atlantic between the Cenomanian and the Maastrichtian with an associated shift from deep-water formation in the South Pacific (Cenomanian) to the South Atlantic and Indian Ocean (Maastrichtian). These changes are associated with a reversal of deep-water fluxes across the Caribbean Seaway between North and South America, which provides a possible explanation for decreasing ε_Nd values throughout the Atlantic during the Late Cretaceous (Donnadieu et al., 2016). That study further suggested that paleogeographic evolution during the Late Cretaceous eliminated conditions necessary for the occurrence of OAEs (Donnadieu et al., 2016).

In this contribution, we present a comparison of Cenomanian and Maastrichtian simulations (as well as a number of sensitivity experiments) similar to the Donnadieu et al. (2016) study but using a recent and higher-resolution earth system model. We evaluate the effect of changes in the depth of major Maastrichtian gateways including the Labrador Seaway, Drake Passage, Caribbean Seaway, and Neotethys Seaway, as well as the effect of decreasing the atmospheric CO₂ concentration. The paper is organized as follows. First, we briefly review the paleogeographic history of major Late Cretaceous gateways to describe the rationale behind prescribed gateway changes. We then explore the evolution of the global ocean circulation in the Late Cretaceous with a particular focus on the changes in intermediate and deep-water currents across the globe. Finally, we compare our simulated ocean circulation with compilations of geochemical data in order to provide an updated picture of the global ocean circulation at the close of the Mesozoic Era.

Figure 1 Bathymetry of the Cenomanian and Maastrichtian configurations and enlargements of regions where bathymetric changes were made in the gateway sensitivity experiments.

2 Paleogeographic considerations

Advances in the knowledge of the geological history of ocean gateways, combined with modeling of the likely effects of those changes, may provide critical arguments in favor of some modes of Late Cretaceous ocean circulation over others. This section summarizes observations on Late Cretaceous paleogeography in critical regions relative to the Cenomanian (∼95 Ma) and early Maastrichtian (∼70 Ma) paleogeographic reconstructions used in our model simulations and sensitivity experiments. These two reconstructions are based on proprietary paleogeographies provided by Getech Plc (Fig. 1), which have been introduced by Lunt et al. (2016) and Tabor et al. (2016).

2.1 Equatorial Atlantic

Rifting between western Africa and eastern Brazil began during the Early Cretaceous (Mascle et al., 1988). Marine waters invaded this narrow corridor from both ends during the Early Aptian, and a shallow marine connection between the Central and South Atlantic oceans existed around 104 Ma (Brownfield and Charpentier, 2006; Ye et al., 2017). The NE–SW motion of the South American plate relative to the African plate was accommodated across transform-related marginal ridges dividing the equatorial Atlantic Ocean into narrow basins during Albian–Cenomanian (Basile et al., 2005; Jones et al., 2007), which restricted seawater exchanges between the Central and South Atlantic oceans and favored euxinic conditions and black shale deposition in these basins (Pletsch et al., 2001; Ye et al., 2017). Deep-water exchange among basins remained limited from the Turonian to the middle Coniacian (Pletsch et al., 2001), but the disappearance of black shales in the equatorial Atlantic during the Campanian suggests the initiation of a reliable supply of oxygenated deep water from the South Atlantic Ocean at this time (Jones et al., 2007), thereby marking the beginning of a fully opened connection between the Central and South Atlantic oceans.

Our Cenomanian and Maastrichtian paleogeographies (Fig. 1) are consistent with this geological history of the equatorial Atlantic Seaway. In our Cenomanian paleogeography, this gateway is restricted to a narrow channel with a maximum depth of ∼2000 m, whereas in the Maastrichtian paleogeography, the Atlantic is opened to full deep-water connection (> 3000 m) between the North and South Atlantic.

2.2 Labrador and East Greenland Seaways

Rifting in the Labrador Sea began during the Early Cretaceous, possibly as early as the Valanginian (Dickie et al., 2011), but the onset of sea-floor spreading took place between the Campanian and the Danian (Roest and Srivastava, 1989; Chalmers and Pulvertaft, 2001). This onset was associated with a deepening of the Labrador Sea as indicated by the presence of agglutinated foraminifera from the Maastrichtian onwards (Kuhnt et al., 1989; Setoyama et al., 2017). East of Greenland, the subsidence of the shallow seas occurred later during the Paleocene (Gernigon et al., 2019).

The proto-Labrador Sea is closed in our Cenomanian paleogeography (Fig. 1). Although evidence suggests that a proto-Labrador Sea potentially existed before the Campanian (Dickie et al., 2011), it would have been restricted to shallow depths with limited influence on interbasinal exchange due to the absence of a northward connection to the Arctic Ocean. The configuration of the proto-Labrador Sea in our Maastrichtian paleogeography (Fig. 1) is in line with the distribution of agglutinated foraminifera (Setoyama et al., 2017), with shallow seas east of Greenland and a deeper proto-Labrador Sea to the south. However, the exact paleodepth of the Maastrichtian Labrador and East Greenland seas is still poorly constrained. We investigate the possibility of the existence of deeper marine channels in the Maastrichtian northern North Atlantic by deepening the Labrador and East Greenland seas to 4000 m. This sensitivity experiment represents an end-member of the deepest paleogeographic configuration of the northern North Atlantic in the Maastrichtian, and we note that a deep East Greenland Sea is not supported by Cretaceous geologic evidence.

2.3 Drake Passage

The history of Drake Passage is intertwined with the evolution of the South America–Antarctic Peninsula–Scotia plate system (Eagles, 2016). The geometrical arrangement of southern South America and the Antarctic Peninsula (AP) has been a matter of debate since the pioneering work of Wegener (1924). Paleomagnetic inclinations measured in rocks from the Fuegian Andes have been shown to be statistically indistinguishable from those of the Antarctic Peninsula for the Late Cretaceous (Poblete et al., 2016), suggesting that the tip of the AP remained close to Tierra del Fuego. In addition, rocks from the Navarino microplate (Fuegian Andes) recorded a 100^∘ counterclockwise rotation over the last 120 Myr, which suggests that the AP and the southern Andes formed a linear and continuous margin during the Early Cretaceous (Poblete et al., 2016). Likewise, a clockwise rotation was found in the apparent polar wander path of the AP coeval with the rotation of the Navarino microplate, thus confirming the oroclinal bending of the Fuegian Andes (Milanese et al., 2019). During the Cenomanian the oroclinal bending was at an early stage, such that the tip of the South American plate was still connected to the AP, with the possible existence of a land bridge allowing terrestrial and fresh water vertebrate taxa interchange (Poblete et al., 2016). The presence of a land bridge for terrestrial exchange does not exclude the possible existence of seawater connections but indicates that any connections would have been restricted to shallow water depths. The oroclinal bending continued during the Late Cretaceous, but the AP and southernmost South America remained close to each other during the Maastrichtian. This geography is supported by paleontological evidence placing the onset of terrestrial faunistic isolation in South America in the Late Paleocene around 58 Ma (Reguero et al., 2014). The final disruption of the AP–Patagonia system occurred during the Early Eocene, but the development of deep-water exchange through the Drake passage only began during the Late Eocene (Scher and Martin, 2006; Lagabrielle et al., 2009). In summary, although the complexity of the South America–Antarctic Peninsula–Scotia plate system's geologic history still hampers comprehensive tectonic reconstructions of the Drake Passage region during the Late Cretaceous, recent evidence indicates that any potential seawater connection would have been restricted to shallow water.

In our Cenomanian paleogeography, the deepest part of the Drake Passage reaches ∼800 m along a narrow corridor, while in the Maastrichtian, only upper ocean water exchange is possible through the Drake Passage as its deepest part reaches ∼400 m. Thus, although our Cenomanian Drake Passage is relatively deep, our Cenomanian and Maastrichtian paleogeographic reconstructions are broadly consistent with paleomagnetic and paleontological data (Fig. 1). However, alternative paleogeographic reconstructions exist, in which the Drake Passage exhibits an even deeper configuration (Sewall et al., 2007; Donnadieu et al., 2016; Niezgodzki et al., 2017). Because the recent study of Donnadieu et al. (2016) documents only minor changes to the global ocean circulation for depths of the Drake Passage lower than 1000 m, we have chosen to prescribe a full deep-ocean connection (4000 m) in our sensitivity experiment in order to maximize the potential impact of the deepening of Drake Passage on ocean circulation, even if these abyssal depths are probably exaggerated for the Maastrichtian (Fig. 1).

2.4 Caribbean Seaway

The Caribbean region has a complex geological evolution, which started during the Jurassic with the dislocation of Pangaea (Pindell and Kennan, 2009). Rifting between North and South America during the Jurassic and Early Cretaceous led to the opening of the proto-Caribbean Seaway. To the west, subduction of the Farallon plate beneath the proto-Caribbean plate during the Early Cretaceous formed an oceanic volcanic arc stretching from the northwestern tip of South America to the southern tip of North America (Pindell and Kennan, 2009). Emplacement of the Caribbean Large Igneous Province (CLIP) starting in the Cenomanian marked a turning point in the history of the Caribbean region. This large (4×10⁶ km³) basaltic oceanic plateau was formed from 94 to 89 Ma (Andjić et al., 2019, and references within) or 95 to 83 Ma (Dürkefälden et al., 2019) from volcanism driven by melting during the initial plume head stage of the Galapagos hotspot. The CLIP was initially located along the southern edge of the North American plate and the northwestern edge of the South American plate, westward of the oceanic island arc (Andjić et al., 2019). Constructed from 8 to 20 km of thick but warm (buoyant) basaltic flows emplaced on the oceanic crust of the Farallon plate, the CLIP prevented subduction of the Caribbean plate (Pindell and Kennan, 2009). During the Cenomanian, the CLIP was located in the Caribbean Seaway, and its buoyancy restricted exchange to shallow water passages (Buchs et al., 2018), including local subaerial emergence, as indicated by volcaniclastic deposits exposed in the Western Cordillera of Colombia (Buchs et al., 2018). The CLIP then progressively moved eastward relative to the North and South American plates during the Late Cretaceous, and new subduction zones were initiated on both the east and west sides of the CLIP, leading to new volcanic oceanic arcs (Pindell and Kennan, 2009). Paleontological evidence also supports restricted water exchange between the Pacific and the Atlantic during the latest Cretaceous (Iturralde-Vinent, 2006; Ortiz-Jaureguizar and Pascual, 2011). Recent research, therefore, suggests that the Caribbean Seaway was relatively shallow across the Late Cretaceous interval.

Our Cenomanian and Maastrichtian paleogeographies are consistent with a shallow Caribbean Seaway. The Cenomanian Caribbean Seaway is deeper than that of the Maastrichtian, which is in reasonable agreement with the progressive formation of the CLIP and its eastward motion between the North and South American plates (Buchs et al., 2018). Although geologic evidence does not support the existence of a deep-water connection between the Pacific and the Atlantic in the Late Cretaceous, alternative paleogeographic reconstructions have been employed, in which the Caribbean Seaway is opened to deep flow (Sewall et al., 2007; Donnadieu et al., 2016). As we did for the Drake Passage, we investigate the consequences of prescribing a full deep-ocean connection through the Caribbean Seaway, by deepening the southern portion of this seaway to 4000 m (Fig. 1).

2.5 Neotethys Seaway

The Neotethyan Ocean exhibits a complex geological history. There is evidence for Late Cretaceous marine exchange between the Central Atlantic Ocean and the Neotethyan Ocean, which mostly occurred through narrow and deep corridors (Stampfli, 2000; Stampfli and Borel, 2002; Nouri et al., 2016). These corridors formed during the final break-up of the Pangaea supercontinent, which led to the opening of the Alpine Tethyan Ocean during the Early Jurassic coeval with the opening of the Central Atlantic Ocean (Stampfli and Borel, 2002). The Alpine Tethyan Ocean began to close in the Early Cretaceous in response to the rotations of the African plate and the Iberian plate (Stampfli and Borel, 2002). During the Late Cretaceous, two deep marine corridors located on both sides of the Anatolides–Taurides permitted water exchanges between the Central Atlantic Ocean and the Neotethyan Ocean (Stampfli and Borel, 2002; Nouri et al., 2016), but it is unclear whether bathymetric sills locally restricted these exchanges to shallow depths (Stampfli and Borel, 2002).

In our paleogeographic reconstructions, the Cenomanian Neotethys allows a deep-water marine connection between the eastern Neotethys and the North Atlantic, whereas the Maastrichtian Neotethyan Ocean does not (Fig. 1). The continued convergence of the African and Eurasian plates throughout the Late Cretaceous (Stampfli and Borel, 2002) can be tentatively used to support the existence of deeper connections in the Cenomanian than the Maastrichtian, but existing uncertainties still preclude any firm conclusions on the absence of deep-water connection through the Neotethyan Ocean in the Maastrichtian. Here, as above, we investigate the consequence of a full deep-ocean connection (4000 m depth) between the Neotethyan Ocean and the North Atlantic.

3 Model description and experimental design

The simulations are performed with the CCSM4 earth system model (ESM) (Gent et al., 2011, and references therein). Our CCSM4 setup is comprised of the POP2 dynamic ocean model, the CAM4 atmosphere model, the CLM4 land surface model, and the CICE4 sea-ice model. The atmosphere and land-surface components run on a finite-volume grid at $\mathrm{1.9}{}^{\circ }\times \mathrm{2.5}{}^{\circ }$ resolution with 26 uneven vertical levels, while the ocean and sea-ice components run on a rotated-pole distorted grid at roughly 1^∘ resolution with 60 vertical levels that vary in thickness with depth.

We perform two baseline simulations of the Cenomanian and early Maastrichtian, which are branched from the 1500-year-long CEN and MAA simulations described in Tabor et al. (2016). The simulations are, respectively, run for 500 and 850 additional years with prescribed vegetation fields adapted from Sewall et al. (2007) rather than the dynamic vegetation model of Tabor et al. (2016) because the latter produced low vegetation density at high latitudes that did not agree well with fossil-based reconstructions (Tabor et al., 2016). As a result, simulated high-latitude land surface temperatures tended to be too cold and seasonally variable. Switching to prescribed vegetation, based on Sewall et al. (2007), helped reduce this temperature bias. Other boundary conditions do not change: the atmospheric CO₂ concentration is set to 1120 ppm (4 times the preindustrial atmospheric levels, PAL =280 ppm) in line with proxy-based reconstructions (Wang et al., 2014). Other greenhouse gas concentrations are set to their preindustrial values. We use a modern Earth orbital configuration, and the total incoming solar irradiance is reduced to appropriate Cenomanian and Maastrichtian values of 1353.9 and 1357.18 W m⁻², respectively, following Gough (1981).

The gateway sensitivity experiments, in which either the Labrador Seaway, Drake Passage, Caribbean Seaway, or Neotethys Seaway is deepened to 4000 m, are branched from the 850-year-long extension of our Maastrichtian simulation. Note that we refer to these bathymetric regions as gateways (or seaways) for simplicity, although they may not be gateways in its truest sense (i.e., a narrow passage connecting two otherwise separated ocean basins). The baseline Maastrichtian case and four sensitivity experiments are each run for another 950 years. We also perform another sensitivity experiment in which atmospheric CO₂ levels are decreased to 560 ppm (2× PAL) in a Maastrichtian configuration because proxy records indicate that the latest Cretaceous was a time of lower CO₂ concentrations than the Cenomanian (e.g., Breecker et al., 2010; Wang et al., 2014; Foster et al., 2017). As the Cenomanian and baseline Maastrichtian simulations, the 2× CO₂ simulation is branched from the 1500-year-long MAA2x simulation described in Tabor et al. (2016) and is run for 1350 years with the Sewall et al. (2007) prescribed vegetation fields and keeping the other boundary conditions identical to that of the baseline Maastrichtian experiment. In total, the Cenomanian and 2× CO₂ Maastrichtian simulations have been run for 2000 and 2850 years, respectively, whereas the baseline Maastrichtian and the gateway sensitivity simulations have been run for 3300 years. Note that the baseline Maastrichtian and 2× CO₂ Maastrichtian simulations are identical to those used in Haynes et al. (2020).

Figure 2 (a) Time series of temperature at the sea surface and in the intermediate (1000 m) and deep (3000 m) ocean, showing that the model has reached quasi-equilibrium at the end of the simulations. The gap between years 850 and 930 in the Maastrichtian simulation (red line) is due to unfortunate loss of data. Only the ends of the sensitivity simulations (Deep Labrador Seaway, Deep Drake Passage, Deep Caribbean Seaway, and Deep Neotethys Seaway) are shown because the full history of the evolution of these simulations was not conserved. Note that the first 1500 years of the simulations, described in Tabor et al. (2016), are omitted on this figure. (b) Time series of the meridional overturning circulation. Note that the maximum overturning intensity is negative because the circulation is counterclockwise.

At the end of integration, the simulations have reached quasi-equilibrium in the deep ocean, as characterized by time series of temperature and meridional overturning circulation (MOC; Fig. 2). A small residual trend exists in the intermediate ocean of the Maastrichtian simulation (1000 m temperatures), which is probably linked to the interval of MOC intensification in this simulation (Fig. 2). This small trend is unlikely to affect the outcomes of this study because the patterns of the ocean circulation do not change during the interval of lower MOC intensity.

Our version of CCSM4 incorporates an ideal age tracer of water masses, a common tool for tracking water mass pathways. Ideal age is an excellent tracer in a fully equilibrated ocean, but for an ocean that is initiated from an unequilibrated state (as in this study), the ideal age tracer is affected by the spinup history and does not track the equilibrium circulation. To use this tracer quantitatively in this study, the simulations would require an additional 2000 years of integration, a computational expense that we could not afford. Alternative techniques, such as Newton–Krylov solvers, exist to estimate the equilibrium values of ocean tracers in an offline procedure (e.g., Bardin et al., 2014; Lindsay, 2017), and such analyses will be the focus of future work. In this paper, we use the ideal age tracer only as a qualitative, complementary diagnostic of deep-water formation regions.

Results presented in the following sections are averaged over the last 100 years of each simulation. We first describe general characteristics of the surface climate, the global overturning circulation, and how ocean temperatures respond to changing paleogeography. Next, we focus on the intermediate and deep circulation and analyze how circulation patterns differ between the Cenomanian and the Maastrichtian. To characterize differences, we track the exchange of water across major oceanic sections by calculating positive and negative water fluxes (Table S1 in the Supplement) for three depth ranges– upper ocean (< 500 m depth), intermediate ocean (500–1500 m), and deep ocean (> 1500 m). Note that we refer to the net exchange across a section as the sum of positive and negative fluxes across the section. Simulated changes in ocean circulation between the Cenomanian and the Maastrichtian and between the Maastrichtian and the sensitivity experiments are then compared to previous modeling studies and geochemical data.

Figure 3 Climate diagnostics for the Cenomanian simulation. (a) Surface ocean (upper 100 m) temperature (^∘C). (b) Surface ocean (upper 100 m) salinity (PSU). (c) Precipitation minus evaporation (mm d⁻¹). (d) Runoff freshwater flux (mSv). (e, f) Late winter maximum mixed layer depth (m).

4 Results

4.1 Cenomanian circulation

4.1.1 Surface climate and global overturning circulation

The global-average annual surface ocean (upper 100 m) temperature of the Cenomanian simulation reaches 26.1 ^∘C. Maximum upper ocean temperatures of more than 34 ^∘C are found in the low-latitudes in the western Pacific Ocean and in the Saharan epicontinental sea in Africa, whereas the eastern Pacific Ocean is much cooler because of wind-driven upwelling (Fig. 3a). Relatively warm (> 10 ^∘C) waters exist in the high latitudes in the Southern Ocean and the North Pacific, though high-latitude coastal and Arctic Ocean waters are colder. Arctic Ocean mean surface ocean temperatures average 2.7 ^∘C. The cold conditions in the Arctic Ocean allow for the formation of winter sea ice (Fig. S1 in the Supplement). The Southern Ocean does not freeze seasonally with the exception of an inlet between the Antarctic and Australian continents (Fig. S1).

The modeled upper ocean salinity generally correlates with patterns of precipitation minus evaporation (PME). The highest open ocean salinities are found in subtropical evaporative areas in the center of major ocean gyres, while lower values are found in the equatorial Neotethyan Ocean and western Pacific and in the high latitudes (Fig. 3b and c). The Arctic Ocean contains low salinity values reflecting the fact that it is a nearly enclosed basin in a region of net freshwater input. In addition, freshwater input from continental rivers affects the spatial distribution of salinity (Fig. 3d), in particular in coastal areas and epicontinental seas. The epicontinental northwestern part of Asia is a region of low salinity due to the isolation of this seaway from the open ocean and of the supply of freshwater from runoff and precipitation. Other low-salinity coastal waters are found at equatorial latitudes in enclosed epicontinental basins on the eastern coast of West Africa (Saharan epicontinental sea) and on the northwestern coast of South America as well as in the isolated high-latitude basin located between Australia and Antarctica. In contrast, high-salinity waters are found in enclosed subtropical basins, such as on the western coast of South America and on the Asian margin of the Neotethyan Ocean as well as in the Gulf of Mexico (Fig. 3b). These high-salinity areas correlate with regions of high temperature, low river freshwater input, and largely negative PME (Fig. 3a–d).

Figure 4 Global meridional overturning circulation (clockwise positive) for each experiment. (a) Cenomanian, (b) Maastrichtian, (c) Deep Labrador Seaway, (d) Deep Drake Passage, (e) Deep Caribbean Seaway, (f) Deep Neotethys Seaway, and (g) 2× CO₂ Maastrichtian.

The Pacific sector of the Southern Ocean is comparatively warmer and more saline than other high-latitude regions. Cooler and fresher waters in the North Pacific are due to the mixing of North Pacific waters with cold, fresh Arctic waters across the Cenomanian Bering Strait. In the Indo-Atlantic sector of the Southern Ocean, seawater salinities are lower than in the Pacific sector due to the large relative freshwater flux from riverine input into a smaller basin (Table S2). The other major reason for this South Pacific anomaly is a temperature- and salt-advection feedback linked to the winter deepening of the mixed-layer depth (MLD) associated with a large area of deep-water formation (Fig. 3e). The same process occurs in the North Atlantic, albeit at a smaller scale in terms of areal extent and of depth reached by sinking waters (Fig. 3f). Predicted global MOC is essentially fed by sinking South Pacific waters, which drive a strong overturning cell in the Southern Hemisphere, with a maximum of ∼18 Sv (sverdrups) around 40^∘ S and 2000 m, and whose lower limb extends to approximately 40^∘ N at depths of ∼4000 m (Fig. 4a). In the Northern Hemisphere, the formation of intermediate waters in the North Atlantic leads to a weak Atlantic Meridional Overturning Circulation (AMOC; Figs. 4a and S2), which reaches up to ∼1500–2000 m around 40^∘ N (Fig. S2).

Figure 5 Intermediate (500–1500 m) water flow across major oceanic sections defined in Table S1 and represented in green for each simulation. Because flows vary across the range of depth depicted, summed flow in both directions across each section is shown with the larger flux in red and the smaller in orange. Thus, the direction of the red arrow gives the direction of the net intermediate flow across a section, and the magnitude of the net flow is given by the difference between the fluxes represented by the red and orange arrows. Abbreviated sections: C (Caribbean), CA (Central Atlantic), D (Drake), EI (East Indian), IA (Indo-Asian), Ind (Indonesian), Med (Mediterranean), SA (South Atlantic), SAf (South African), SC (South China), Tet (Tethys), and WI (West Indian).

4.1.2 Intermediate (500–1500 m) circulation

The intermediate ocean circulation is fed primarily by two sources: upwelling of deep waters and sinking of upper ocean waters to intermediate depths (Fig. 3e and f). North Atlantic intermediate waters are composed of upper waters that sink in the North Atlantic, upwelled deep waters from the Neotethyan Ocean that were advected across the Mediterranean (Table S1, Mediterranean section, and Fig. 5a), and weaker inputs of intermediate waters from the Pacific and Central Atlantic (Table S1, Caribbean and Central Atlantic sections, and Fig. 5a). More than 90 % of the intermediate waters advected out of the North Atlantic flow westward across the Caribbean Seaway (Table S1, Caribbean and Central Atlantic sections), while the remaining fraction flows southward through the South Atlantic to the Southern Ocean where it is joined by weak transport of Pacific intermediate waters across the Drake Passage (Table S1, South Atlantic and Drake sections, and Fig. 5a). The South Atlantic intermediate waters are then advected eastward trough the Indian Ocean and then northeastward to the eastern Neotethyan Ocean (Figs. S3 and S4). These waters eventually flow into the Pacific by joining an eastern Neotethyan Ocean recirculation of Pacific intermediate waters, forming a narrow, intense eastward current that follows the Australian coast (Fig. S4). In contrast, intermediate waters circulating northwestward toward the western Neotethyan Ocean originate exclusively from Pacific intermediate waters. The Pacific intermediate water system is essentially comprised of a mixture of North Atlantic intermediate waters that flow westward through the Caribbean Seaway, upwelled Pacific deep waters, and recirculated Indian Ocean intermediate waters mentioned above.

Figure 6 Deep (> 1500 m) water flow across major oceanic sections defined in Table S1 and represented in green for each simulation. As in Fig. 5, because flows vary across the range of depth depicted, summed flow in both directions across each section is shown with the larger flux in red and the smaller in orange. Thus, the direction of the red arrow gives the direction of the net intermediate flow across a section and the magnitude of the net flow is given by the difference between the fluxes represented by the red and orange arrows. Abbreviated sections: C (Caribbean), CA (Central Atlantic), D (Drake), EI (East Indian), IA (Indo-Asian), Ind (Indonesian), Med (Mediterranean), SA (South Atlantic), SAf (South African), SC (South China), Tet (Tethys), and WI (West Indian).

Figure 7 Cenomanian deep circulation (3000 m) in (a) the eastern Neotethyan Ocean and (b) the southwestern Pacific Ocean. Orange arrows represent major deep current systems in the eastern Neotethyan Ocean and southwestern Pacific Ocean. Purple contours represent regions of deep waters formation (contour 500 m). Section A–B defines the Indonesian section of Table S1. (c) Fluxes of water across the Indonesian section over the whole water column.

4.1.3 Deep (> 1500 m) circulation

The southwestern Pacific is the source region for deep waters in our Cenomanian simulation (Fig. 3e). These sinking waters either fill the deep eastern Pacific basin or are advected westward across the Indonesian section, following a strong coastal current around Australia (Figs. 6a and 7). Deep waters crossing the Indonesian section following this westward current mostly recirculate back to the Pacific and mix with the eastern Pacific deep waters to fill the North Pacific basin. Less than 10 % of the westward-flowing deep waters that have crossed the Indonesian section are advected southward across the East Indian section to the Southern Ocean (Table S1, Indonesian and East Indian sections). Deep waters exported northward toward to western Neotethyan Ocean mostly come from a deep intermediate westward current that follows the southern tip of Asia between ∼800 and 2400 m (Fig. 7c and Table S1, Indo-Asian and Tethys sections). In the Southern Ocean, deep waters are advected to the South Atlantic, but regions of shallow bathymetry (e.g., the Kerguelen Plateau) largely restrict deep-water flow, and these waters ultimately well up to shallower depths (Fig. 6a). The fate of deep waters flowing northward from the Neotethyan basin is similar. These waters are advected across the western Neotethyan Ocean to the North Atlantic, where they are upwelled to shallower depths because the Caribbean Seaway is closed to deep flow (Fig. 6a). An examination of the zonally averaged ideal age values in the Atlantic basin reveals that the deepest-sinking waters in the North Atlantic winter MLD regions reach the deep ocean (Figs. S5 and 3f). These waters are mostly restricted to the North Atlantic; indeed, only a tiny fraction of North Atlantic deep waters is advected southward into the Central Atlantic (Table S1).

In summary, the bathymetric restrictions in the Cenomanian Atlantic, western Neotethys and Southern Ocean largely confine deep-water circulation to the Pacific and eastern Neotethyan Ocean. In contrast, a vigorous intermediate circulation marked by a strong circum-equatorial global current exists, although the restricted Central and South Atlantic basins remain mostly stagnant.

Figure 8 Climate diagnostics for the Maastrichtian simulation. (a) Surface ocean (upper 100 m) temperature (^∘C). (b) Surface ocean (upper 100 m) salinity (PSU). (c) Precipitation minus evaporation (mm d⁻¹). (d) Runoff freshwater flux (mSv). (e, f) Late winter maximal mixed layer depth (m).

Figure 9 Climate diagnostics for the Maastrichtian simulation relative to the Cenomanian simulation. (a) Surface ocean (upper 100 m) temperature difference (^∘C). (b) Surface ocean (upper 100 m) salinity difference (^∘C). (c, d) Cenomanian and Maastrichtian barotropic streamfunction (Sv).

4.2 Maastrichtian circulation

4.2.1 Evolution of surface climate and global overturning circulation

The combined changes in paleogeography and solar constant from the Cenomanian to the Maastrichtian lead to a global SST warming of only ∼0.1 ^∘C, suggesting that changes in paleogeography may cause cooling that compensates for the increasing solar constant (Lunt et al., 2016). Though the global temperature change is minimal, there are substantial regional temperature and salinity changes at the surface in the Maastrichtian compared to the Cenomanian (Tabor et al., 2016). Maastrichtian North Pacific surface ocean waters warm significantly because of the closure of the Arctic connection (Figs. 8a and 9a). As a result, the Arctic Ocean becomes more enclosed, cools, and, as a region of net freshwater input, freshens (Figs. 8a–d and 9a–b). The reduction in the intensity of the circum-equatorial global current (Table S1 and Fig. 9c–d) in the Maastrichtian reduces coastal upwellings of deeper and colder waters on the northern coast of Africa and South America, leading to surface warming of up to a few degrees. The eastern equatorial Pacific warms because of a weaker Walker circulation, which reduces the east–west ocean temperature gradient (Poulsen et al., 1998; Tabor et al., 2016). The PME in the eastern equatorial Pacific increases in the Maastrichtian and leads to lower salinity (Figs. 8b–c and 9b). The opening of the South Atlantic Ocean and Southern Ocean during the Late Cretaceous created a wider basin, which allows for a large subpolar gyre to form (Fig. 9d). This gyre reduces the advection of warm and saline subtropical waters into the Southern Ocean along the eastern coast of Africa (Fig. 9c and d), and the Southern Ocean cools and freshens as a result (Fig. 9a and b). In addition, the Ekman pumping associated with the subpolar gyre leads to upwelling of deeper and colder waters to the surface, which contributes to cooling of the South Atlantic and Indian oceans. In the northern part of the eastern Neotethyan Ocean, the salinity increase (Fig. 9b) is due to changes in the patterns of surface currents, which limits the northward advection of fresher Neotethyan equatorial waters in the Maastrichtian. Finally, cooling in the North Atlantic and warming in the Pacific sector of the Southern Ocean are related to changes in the MOC (Fig. 4b). In contrast to the Cenomanian, the Maastrichtian North Atlantic does not exhibit deep intermediate water formation (Fig. 8e and f). This elimination of proto-AMOC weakens the advection of warm and saline subtropical waters into the North Atlantic, leading to surface cooling. Conversely, the intensification of South Pacific deep-water formation drives a more expansive global MOC (Fig. 4b) and is associated with surface warming (Fig. 9a) via reinforcement of the temperature- and salt-advection feedback.

Figure 10 Zonally averaged temperature difference (^∘C) between the Maastrichtian and the Cenomanian simulations. (a) Global average. (b, c, d) Pacific, Atlantic, and Indo-Neotethyan average, based on basins defined in Fig. S3.

4.2.2 Temperature changes in the intermediate and deep ocean

The Cenomanian to Maastrichtian paleogeographic evolution, in particular the widening of the Atlantic Ocean, the northward migration of the Indian and Australian subcontinents, and the varying configuration of major gateways, results in a complete reorganization of intermediate and deep ocean circulation (Table S1 and Figs. 5a–b and 6a–b). This reorganization leads to significant changes in temperatures in the ocean interior in the Maastrichtian relative to the Cenomanian (Figs. 10 and S6). The global temperature change essentially reflects the Pacific signal because of the size of the Pacific basin in both stages (Fig. S3). In the South Pacific Ocean, increased ventilation in the Maastrichtian explains most of the warming signal (Figs. 10b and S6). In the North Pacific, Maastrichtian intermediate water cooling is attributed to restriction to shallow water depths (< 500 m) of flow through the Caribbean Passage, hampering westward advection of North Atlantic waters below the uppermost ocean layers. It is important to note that this restriction is only significant because, in the Cenomanian, North Atlantic intermediate waters are warmer than North Pacific intermediate waters due to deep-water formation occurring in the North Atlantic. In the Maastrichtian, due to the absence of deep-water formation in the North Atlantic, intermediate waters are colder because of reduced ventilation (Figs. 10c and S6) and the geometry and bathymetry of the Caribbean Seaway and of the western Neotethyan Ocean, which isolates the basin from intermediate and deep-water exchange with the Pacific and Neotethyan oceans.

The northward displacement of India and the widening of the Atlantic in the Maastrichtian paleogeography reduce the isolation of the deep South Atlantic, and this basin is invaded by deep flow from the Pacific via the Indian Ocean (Table S1, East Indian and South African sections, and Fig. 6b) leading to lower temperatures (Figs. 10c and S6). Finally, the Neotethyan basin is mostly warmer in the Maastrichtian than it is during the Cenomanian (Fig. 10d). The small deep ocean warming is explained by advection of warmer deep waters formed in the South Pacific. The larger upper intermediate ocean (centered on ∼500 m depth) warming is explained by differences in the configuration of the western Neotethyan Ocean. In the Cenomanian simulation, Neotethyan upper intermediate waters, formed in the late winter when the MLD deepens (Fig. 3f), are advected toward the North Atlantic because the western Neotethyan Ocean is open to intermediate and deep waters (Fig. S7a). The closure of the western Neotethyan Ocean to intermediate and deep waters in the Maastrichtian simulation hampers this advection, and flow of these waters shifts southward to the eastern Neotethyan Ocean (Fig. S7b). These sinking upper intermediate waters carry a higher temperature and salinity signal into the eastern Neotethyan Ocean, which can be followed on transects across the Neotethyan Ocean (Fig. S7c–f), and are responsible for the warming of this part of the basin in the Maastrichtian.

4.2.3 Evolution of the intermediate (500–1500 m) circulation

With the restriction of intermediate and deep flow through the Caribbean Seaway and the western Neotethyan Ocean, the sources of intermediate waters in the North Atlantic Ocean are deep waters advected from the South Atlantic that are upwelled in the North Atlantic (Fig. 5b) and winter downwelling of upper ocean waters in the northern part of the basin (Fig. 8f). North Atlantic intermediate waters return to the Pacific via the South Atlantic and the Indian Ocean (Table S1), following a strong eastward coastal current around the northern tip of Australia (Fig. S8) similar to that existing in the Cenomanian simulation (Fig. S4). In the eastern Neotethyan Ocean, intermediate waters are primarily composed of intermediate Pacific waters that flow westward across the Indonesian section between 0 and 10^∘ S (Fig. S8), of eastern Neotethyan Ocean deep waters that are upwelled to shallower depths (Table S1, Indo-Asian section, and Fig. 6), and of winter upper ocean waters that sink in the northern part of the eastern Neotethys (Fig. S7b). These eastern Neotethyan Ocean intermediate waters flow eastward into the Pacific following a southward current along the eastern Indian margin and mostly join the strong eastward current circulating around Australia (Fig. S8).

Figure 11 Maastrichtian deep circulation (3000 m) in (a) the eastern Neotethyan and Indian oceans and (b) the southwestern Pacific Ocean. Orange arrows represent major deep current systems in the eastern Neotethyan and Indian oceans and southwestern Pacific Ocean. Purple contours represent regions of deep waters formation (contour 500 m). Section A–B defines the Indonesian section of Table S1. (c) Fluxes of water across the Indonesian section over the whole water column.

4.2.4 Evolution of the deep (> 1500 m) circulation

In the Maastrichtian, as in the Cenomanian, deep waters are formed in the South Pacific, mostly in the western part of the basin, and flow northwestward along the Australian coast (Fig. 11). Along the northern continental slope of the Australian margin, deep waters either cross the Indonesian section eastward or recirculate to fill the Pacific basin (Table S1 and Fig. 11). As in the Cenomanian, deep waters advected across the Indonesian section then either fill the Indian Ocean (Table S1, East Indian section) or journey northward to recirculate toward the Pacific Ocean or the eastern Neotethyan Ocean (Table S1, Indo-Asian section and Fig. 11). Because the connections through the western Neotethyan Ocean are restricted to shallow flow in the Maastrichtian, there is no deep flow across the western Neotethys (Table S1, Tethys and Mediterranean sections, and Fig. 6b). In contrast, the opening of the South Atlantic and Southern Ocean allows stronger deep-water flow from the Indian Ocean into the South Atlantic (Table S1, East Indian, West Indian, and South African sections, and Fig. 6b), which is then advected northward to the North Atlantic (Table S1, South and Central Atlantic sections) and progressively upwelled to shallower depths.

In the Maastrichtian simulation, the net deep circulation appears to flow in the opposite direction of the intermediate circulation (Figs. 5b and 6b). The Maastrichtian circulation is also characterized by more intense meridional exchanges (compare Cenomanian and Maastrichtian meridional sections in Table S1, for instance the East Indian, South Atlantic, and Central Atlantic sections), whereas the Cenomanian circulation is dominated by zonal flow (Table S1, for instance the Indonesian, Tethys, and Caribbean sections).

4.3 Sensitivity of the Maastrichtian circulation to ocean gateways and atmospheric CO₂

As shown above, changes in paleogeography between the Cenomanian and Maastrichtian lead to substantial changes in simulated intermediate and deep ocean circulation. In this section, we analyze the influence of specific gateways and lower atmospheric CO₂ levels on Maastrichtian ocean circulation.

Figure 12 Surface, intermediate, and deep ocean temperature difference (^∘C) between the sensitivity experiments and the Maastrichtian simulation.

4.3.1 Deepening of the Labrador Seaway

Temperature changes in the ocean

Deepening the Labrador Seaway only marginally impacts the global ocean circulation. In this experiment, as in the baseline Maastrichtian configuration, deep-water formation takes place in the South Pacific and mostly in the western part of that basin. The maximum winter MLD in both hemispheres is only weakly different from that of the Maastrichtian (Fig. S9a), and the resulting MOC is nearly identical in structure and intensity (Fig. 4c). In the northern North Atlantic and western Neotethyan Ocean, the slight deepening of the maximum winter mixed layers (Fig. S9a) is associated with surface ocean warming, whereas the surface ocean cools south of Greenland (Fig. 12a). There are only minor temperature changes in other oceanic basins or in the ocean interior (Figs. 12a and S10).

The pattern of upper ocean temperature change is linked to the altered bathymetry of the Deep Labrador Seaway experiment, leading to substantial reorganization of upper ocean currents in the northern North Atlantic (Fig. S11). In the Maastrichtian simulation, waters originating from the North Atlantic subtropical latitudes are largely confined south of Greenland because the shallow bathymetry of the seas bathing the east of Greenland and modern Europe (Fig. S11a). An intense southward flow originating from higher Arctic latitudes exist along the eastern margin of Greenland. This flow then circulates southeastward around the southern edge of the Eurasian continent toward the western Neotethyan Ocean. In the Deep Labrador Seaway experiment, the deepening of the seas south and east of Greenland breaks the confinement of North Atlantic subtropical waters south of Greenland, which are instead advected eastward toward the western Neotethyan Ocean along the southern margin of the Eurasian continental landmass (Fig. S11b). This eastward current also blocks the southern penetration of the east Greenland current originating from Arctic latitudes, the intensity of which is also reduced. In summary, warm subtropical waters flow eastward in the Deep Labrador Seaway experiment rather than being confined south of Greenland, which cools the upper ocean there and warms the western part of Europe. In the region east of Greenland, the decreased supply of cold high-latitude waters leads to warming (Fig. 12a).

Intermediate and deep circulation changes

There are no changes in the direction of intermediate and deep-water transports across major oceanic sections in the Deep Labrador Seaway experiment relative to the Maastrichtian simulation (Figs. 5b–c and 6b–c). The water fluxes are generally slightly higher, which is probably linked to the deepening of the North Atlantic and western Neotethyan oceans winter MLD and associated increase in the vigor of ocean circulation (Fig. S12).

4.3.2 Deepening of the Drake Passage

Temperature changes in the ocean

Deepening of the Drake Passage has a more significant effect on global ocean circulation than the deepening of the Labrador Seaway. Although deep-water formation still occurs in the South Pacific, the intensity of the MOC decreases (Fig. 4d) because deep-water formation is greatly reduced in the South Pacific, in particular along the eastern edge of Zealandia (Fig. S9b). At the latitudes of the Drake Passage, the MLD increases across the whole South Pacific (Fig. S9b) because of the establishment of a deep-water connection through the Drake Passage, which increases the intensity of the eastward current in the South Pacific. The reduction in the intensity of deep-water formation drives upper ocean temperature cooling in the South Pacific, which is partly carried, albeit weakly, at depth to the Atlantic through the Drake Passage (Figs. 12b and S13). The Atlantic is thus better ventilated because the deep Drake Passage connection allows newly formed, young deep waters to invade the Atlantic (Table S1). In contrast, the North Pacific and Neotethyan oceans are less well ventilated because of lower rates of deep-water formation and a lower advection of deep waters across the Indonesian section (Table S1), associated with a small warming.

Intermediate circulation changes

The intermediate circulation with an open Drake Passage undergoes only a few changes relative to the Maastrichtian. An eastward current develops across the Drake Passage and joins the southward flow from the Atlantic Ocean. This increase in the net supply of intermediate waters in the Southern Ocean (Table S1, Drake, South Atlantic, and South African sections) drives a reversal of the intermediate circulation west of India (Table S1, West Indian section, and Fig. 5b and d). This northward water flux enhances the intensity of the intermediate circulation in the eastern Neotethyan Ocean (Table S1, Indo-Asian section), but the structure of the circulation does not change (Figs. S8 and S14). The Pacific intermediate circulation is also similar in the Drake Passage experiment as it is in the Maastrichtian simulation.

Deep circulation changes

The deep circulation in the equatorial eastern Neotethyan Ocean and at the Neotethyan-Pacific boundary does not change (Fig. S14), but opening the Drake Passage to deep circulation significantly reduces the flux of deep water flowing westward across the Indonesian section and into the Indian sector of the Southern Ocean (Table S1 and Fig. 6d). This change is balanced by eastward flow across the Drake Passage, which becomes the dominant source of deep waters in the Atlantic sector of the Southern Ocean. In the Indian and Neotethyan oceans, most of the water flow directions are similar to the Maastrichtian simulation except west of India where the net southward deep-water flow stops. In contrast to the Maastrichtian simulation, with deepening of the Drake Passage, deep waters in the South and North Atlantic mostly originate from Pacific waters flowing eastward through the Drake Passage rather than waters from the Indian Ocean.

4.3.3 Deepening of the Caribbean Seaway

Temperature change in the ocean

Similar to the deepening of the Drake Passage, the opening of the Caribbean Seaway to deep flow causes profound restructuring of the global ocean circulation. Deep-water formation continues to take place in the South Pacific with a reduction in the depth of the winter mixed-layer east of Zealandia relative to the Maastrichtian simulation (Fig. S9c). Consequently, the global MOC is slightly weaker between 2000 and 3000 m (Fig. 4e). The deepening of the Caribbean Seaway leads to cooling of the Atlantic intermediate and deep waters and only minor temperature changes in the Pacific, Indian, and Neotethyan oceans relative to the Maastrichtian, whereas it leads to limited and spatially heterogeneous upper ocean temperature changes (Figs. 12c and S15). As in the Deep Drake Passage experiment relative to the Maastrichtian, the Atlantic Ocean is better ventilated in the Deep Caribbean Seaway experiment than in the Maastrichtian simulation, although intermediate and deep waters invade the Atlantic from the north of the basin rather than from the south.

Intermediate circulation changes

As in the Deep Drake Passage experiment, deepening the Caribbean Seaway does not cause major changes to the modeled global intermediate circulation compared to the Maastrichtian simulation. Changes include the development of weak exchanges of similar magnitude between the Atlantic and the Pacific across the Caribbean Seaway as well as the reversal of the intermediate flow across the West Indian section (Table S1 and Fig. 5e). However, the fluxes of water transported by these altered flows are small, and the overall structure of the intermediate circulation in the Deep Caribbean Seaway remains similar to that of the Maastrichtian (Table S1 and Fig. 5e).

Deep circulation changes

The most salient consequence of the deepening of the Caribbean Seaway on the deep circulation is the reversal of the water fluxes in the Atlantic, from a net northward-dominated flow in the Maastrichtian simulation to a southward-dominated flow in the Deep Caribbean Seaway experiment (Fig. 6b and e) due to the invasion of Pacific deep waters into the Atlantic. In the Southern Ocean, the net transport of water shifts from westward-dominated transport to eastward-dominated transport across the South African section (Table S1 and Fig. 6b and e). As in other Maastrichtian simulations, deep waters formed in the South Pacific flow across the Indonesian section and are either advected into the Indian sector of the Southern Ocean or recirculated to the Pacific (Figs. 11 and S14). However a stronger eastward deep-water flow exists at the southern tip of the Asian continent because of the entrainment created by the opening of the Caribbean Seaway to deep circulation (Table S1, Figs. 6e, and S14). This strong current and the reversal of the net transport of deep waters between the Atlantic and Indian sectors of the Southern Ocean induce a reversal of the deep flow west of India (Table S1, West Indian section, and Fig. 6e). The Southern Ocean is filled with a combination of westward-flowing Indian Ocean deep waters and southward-flowing Atlantic deep waters, which originate from the Pacific and have been advected through the Caribbean Seaway.

4.3.4 Deepening of the Neotethys Seaway

Temperature change in the ocean

In the Maastrichtian and sensitivity simulations described so far, the Neotethys Seaway is shallow and inhibits intermediate and deep ocean circulation (Fig. 1). The deepening of the Neotethys Seaway causes a significant reorganization of the circulation. As in the Deep Drake Passage and Deep Caribbean Seaway simulations, deep-water formation occurs in the South Pacific, although the maximum late winter MLD is reduced relative to the Maastrichtian simulation (Fig. S9d), leading to a slight slowdown of the global MOC (Fig. 4f). Changes in ocean temperatures are minor except in the North Atlantic and Neotethyan oceans at intermediate depth (Figs. 12d and S16). At these depths, the eastern Neotethyan Ocean cools slightly, and the western Neotethyan and North Atlantic warm slightly (Figs. 12d and S16). These changes are due to the opening of intermediate and deep connections between the North Atlantic and Neotethyan oceans. The warmer and saltier sinking winter upper intermediate waters (∼500 m depth) in the eastern Neotethyan Ocean (Fig. S9d) are advected toward the North Atlantic rather than the Indian Ocean (Fig. S17), which leads to the observed intermediate temperature signal. It is noteworthy that this reorganization of water currents caused by the deepening of the Neotethys Seaway is opposite the reorganization caused by the restriction of the Neotethys Seaway that occurs between the Cenomanian and the Maastrichtian (Figs. S7 and S17).

Intermediate circulation changes

In the Deep Neotethys Seaway experiment the directions of the net intermediate transports of water across oceanic sections are also similar to that of the Maastrichtian (Table S1 and Fig. 5b and f). The deep western Neotethyan Ocean provides an outlet for North Atlantic intermediate waters across the Mediterranean section, which increases the intermediate water fluxes out of the North Atlantic (Fig. 5f). However, part of these eastward-flowing intermediate waters recirculate to the North Atlantic, both in the uppermost intermediate ocean (∼500 m), where they join the westward-flowing waters that have sunk in winter in the eastern Neotethyan Ocean (Fig. S17), and in the deeper intermediate ocean (Fig. S18). As a consequence, the net intermediate water transport across the Mediterranean section only slightly increases from 0.2 Sv in the Maastrichtian simulation to 0.5 Sv in the Deep Neotethys Seaway experiment (Fig. 5f). The invasion of the Neotethyan Ocean with North Atlantic intermediate waters also reduces the inflow of Pacific intermediate waters in the eastern Neotethyan Ocean (Table S1, Tethys and Indo-Asian sections, and Figs. 5f and S18), which leads to the reversal of the intermediate flow across the eastern West Indian section (Table S1 and Fig. 5f). Other net intermediate transports remain in the same direction as in the Maastrichtian simulation.

Deep circulation changes

The main circulation difference caused by the deepening of the Neotethys Seaway is a reversal of the deep-water flow direction in the Atlantic basin from northward to southward (Figs. 6b and f). In the equatorial Neotethyan Ocean and Neotethyan-Pacific boundary, deep water circulation is similar to that in the Maastrichtian simulation (Fig. S14); however, the deep eastward Pacific return flow is reduced (Fig. 6f and Table S1, Indonesian section). This change is because the deepening of the Neotethys Seaway opens a deep-water pathway for westward-flowing deep waters formed in the South Pacific. These South Pacific deep waters divide between a southwestward component, which flows into the Indian sector of the Southern Ocean, and a northwestward component, which flows into the Neotethyan Ocean (Fig. 6f). The northwestward deep-water flow across the Neotethyan Ocean induces a reversal of the deep circulation west of India, from a southward-dominated flow in the Maastrichtian to a northward-dominated flow in the Deep Neotethys Seaway experiment. The Neotethyan deep waters then flow into the Atlantic sector of the Southern Ocean via the North Atlantic, which explains the reversal of deep-water flow in this basin. The Southern Ocean is bathed by a combination of deep waters coming from the southern Indian Ocean route and from the Neotethyan-Atlantic route (Fig. 6f).

4.3.5 Decreasing atmospheric CO₂ concentration

Temperature changes in the ocean

Reducing the atmospheric CO₂ concentration only marginally impacts the simulated ocean circulation, even though a ∼2.5–3 ^∘C cooling is observed both at the surface and in the ocean interior (Figs. 12e and S19). Deep-water formation occurs in the South Pacific as in the Maastrichtian and gateway sensitivity experiments (Fig. S9e). Maximum late winter MLD increases in the western part of the South Pacific and decreases in the eastern part relative to the Maastrichtian simulation. The global MOC slightly intensifies (Fig. 4g) because the MLD increase occurs in the region where deepest waters are formed (Fig. 8e).

Intermediate and deep circulation changes

As in the Deep Labrador Seaway experiment, there are no changes in the direction of the intermediate and deep circulation in the 2× CO₂ experiment compared to the Maastrichtian (Figs. 5b and g and 6b and g). The more intense 2× CO₂ MOC drives enhancement of the intermediate and deep fluxes (Table S1) across most oceanic sections, but the absence of significant changes in the water mass pathways indicates that the simulated cooling is exclusively due to the radiative effect of the lower atmospheric CO₂ concentration.

5 Discussion

With the exception of the Deep Labrador Seaway and the 2× CO₂ experiments, each gateway change profoundly alters Maastrichtian deep ocean water mass pathways. The deepening of the Drake Passage and Caribbean and Neotethys Seaways opens barriers to deep circulation, leading to changes in the intensity of circulation and pathways of deep-water flow. At intermediate depths, gateway changes affect the origin and intensity of intermediate circulation but have a lesser effect on the flow pathway within and between basins.

5.1 Comparison to previous model results

5.1.1 Late Cretaceous changes in ocean circulation

To our knowledge only Donnadieu et al. (2016) has investigated changes in ocean circulation from the beginning to the end of the Late Cretaceous. That study uses the FOAM model (Jacob, 1997) to conduct simulations of the Cenomanian/Turonian and Maastrichtian using paleogeographies from Sewall et al. (2007). Donnadieu et al. (2016) (hereafter D16) report that the deep ocean circulation in FOAM is highly sensitive to Late Cretaceous paleogeographic evolution and that these paleogeographic changes are responsible for a shift in the sources of Atlantic deep waters and a reversal of the Atlantic deep-water flow, which provide an explanation for the observed decrease in ε_Nd in the Atlantic and Indian Ocean during the Late Cretaceous. Our simulations differ substantially from those of D16 in the paleogeography employed, in particular the configuration of ocean gateways, and in the locations of deep-water formation, which critically affects the simulated pathways of intermediate and deep water masses.

The baseline Cenomanian simulation of D16 shows deep-water formation in the North and South Pacific as well as the South Atlantic. North Atlantic deep waters are sourced from the Pacific and enter the Atlantic through a relatively deep Caribbean Seaway (2000–2500 m), whereas deep waters formed in the South Atlantic are mostly advected toward the eastern Neotethyan Ocean (Figs. 2a and 3c in D16). In our Cenomanian simulation, in which the Caribbean Seaway is closed to deep flow, North and South Atlantic deep waters originate from the Pacific via Neotethyan and Indian routes, respectively (Fig. 6a).

The baseline Maastrichtian simulation of D16 exhibits a shift in deep-water formation from the South Pacific to the southern Indian Ocean, while deep-water formation in the North Pacific and South Atlantic persists. In those simulations, enhanced South Atlantic deep-water formation drives enhanced northward export of deep waters into the North Atlantic, and these deep waters are advected into the Pacific through a deep Caribbean Seaway (Figs. 2b, 3b and d in D16). In our baseline Maastrichtian simulation, the South and North Atlantic are ventilated by deep waters forming in the South Pacific and flowing westward along a pathway through the Indian Ocean, but the shallow Caribbean and Neotethys Seaways confine deep water in the North Atlantic. Interestingly, our Deep Caribbean Seaway experiment, in which the configuration of the Caribbean Seaway is closer to that of the Maastrichtian simulation of D16, predicts a Pacific to Atlantic flow of deep waters across the Caribbean Seaway (Fig. 6e), whereas D16 experiment predicts the opposite. Contrasts in these model results are directly linked to the different areas of deep-water formation in the Southern Ocean predicted by the two models.

The substantial differences between CCSM4 and FOAM and in the details of the simulations make it difficult to unambiguously explain the substantial changes in the source and circulation of deep waters. In comparison to FOAM, CCSM4 is more complex and has higher spatial resolution. In addition, FOAM and CCSM4 simulations differ in the details of the paleogeographies and initial conditions, which hamper explicit examination of why the two models do not form deep waters in the same locations. However, we speculate that freshwater supply via continental runoff is one mechanism that might lead to these different locations of deep-water formation. In both our Cenomanian and Maastrichtian simulations, the South Pacific is a region of low runoff supply relative to the other sectors of the Southern Ocean (Table S1, Figs. 3d and 8d, and Fig. S20a–b). In addition, the higher elevation and more extensive meridional span of the Rocky Mountains in our reconstructions (Fig. S20c–d) compared to the Sewall et al. (2007) paleogeography used by D16 (Figs. 4 and 5 of Sewall et al., 2007) blocks the advection of moisture across North America (e.g., Maffre et al., 2018), which contributes to decreased surface salinity and prevents deep-water formation in the North Pacific. Finally the lower resolution of FOAM in the atmosphere (7.5^∘ longitude by 4.5^∘ latitude) smooths the Rocky Mountains even more. As a consequence, the moisture flux out of the North Pacific driven by Northern Hemisphere westerlies may be enhanced in D16, leading to increased North Pacific surface salinity and more favorable conditions for deep-water formation.

5.1.2 Sensitivity of ocean circulation to atmospheric CO₂ levels

Ocean circulation is mostly insensitive to reducing the atmospheric CO₂ concentrations in our Maastrichtian configuration. The intermediate and deep water mass pathways are identical, although the intensity of the water fluxes across major oceanic gateways is slightly enhanced in the 2× CO₂ simulation (Haynes et al., 2020). This insensitivity of Late Cretaceous ocean circulation to CO₂ levels is consistent with the results of Donnadieu et al. (2016), which shows that Late Cretaceous simulations performed at 2×, 4×, and 8× CO₂ PAL predict similar areas of deep-water formation. In contrast, Farnsworth et al. (2019) recently reported that reducing atmospheric CO₂ levels from 4× to 2× in a Maastrichtian configuration in the HadCM3BL-M2.1aD earth system model led to a shift in deep-water formation area from the South Pacific Ocean to the South Atlantic and Indian oceans. This high sensitivity to CO₂ only occurs in the Maastrichtian simulation among all the 12 Cretaceous simulations (one per Cretaceous stage) performed by Farnsworth et al. (2019) and occurs again only once (in the Selandian stage, ∼60.6 Ma) among all their seven Paleogene simulations. In the other simulations, both the 2× and 4× CO₂ simulations predict similar areas of deep-water formation. The temporal proximity of the Maastrichtian and Selandian stages led Farnsworth et al. (2019) to suggest that the time period close to the Cretaceous/Paleogene boundary might be particularly sensitive to atmospheric CO₂, but it is not clear in this case why their simulation of the Danian stage (∼63.9 Ma) does not exhibit a similar behavior. As Farnsworth et al. (2019) do not provide a detailed analysis of ocean circulation changes in the Maastrichtian and Selandian stages relative to the others, we can only speculate that these changes might be partly caused by high-latitude smoothing, which is performed on the simulations to ensure model stability and which varies between stages (Lunt et al., 2016; Farnsworth et al., 2019).

More generally, the impact of atmospheric CO₂ levels on ocean circulation has been shown to significantly vary in past greenhouse climate modeling work (Poulsen et al., 2001; Lunt et al., 2010; Poulsen and Zhou, 2013; Donnadieu et al., 2016; Hutchinson et al., 2018; Farnsworth et al., 2019; Zhu et al., 2020). The causes for this large spread in results may be diverse and are difficult to isolate, but we hypothesize that the model climate sensitivity to CO₂ and the range of atmospheric CO₂ levels investigated could explain such variability. Winguth et al. (2010) report results of Paleocene–Eocene Thermal Maximum (PETM) simulations using the CCSM3 fully coupled model (with the CAM3 atmospheric model) and show that ocean circulation and deep-water formation areas remain similar regardless of CO₂, although the intensity of overturning decreases with increasing CO₂. More recently, Zhu et al. (2020) report results of PETM simulations performed at 1×, 3×, 6×, and 9× CO₂ PAL using the CESM1.2 ESM (with the CAM5 atmospheric model) and document a shift in deep and intermediate water formation areas between 1× and 3× CO₂ and complete cessation of deep-water formation at 6× and 9× CO₂. The climate sensitivity of CESM1.2 has been shown to be greater than that of CCSM3 and, contrary to CCSM3, to increase with background CO₂ levels (Zhu et al., 2019). Earth system models with high climate sensitivity to CO₂ may demonstrate a higher sensitivity of ocean circulation to CO₂ because the climate state in which the radiative forcing of CO₂ leads to a warming sufficient to stop deep-water formation can be expected to occur for smaller changes in atmospheric CO₂ levels.

5.2 Evolution of intermediate and deep-water circulation during the Late Cretaceous

5.2.1 Neodymium isotope compilation

A compilation of Cenomanian and Maastrichtian ε_Nd values is shown in Fig. 13 and Tables S3 and S4 (modified from Moiroud et al., 2016). The ε_Nd values at each site are averaged between 100 and 90 Ma for the Cenomanian and between 75 and 65 Ma for the Maastrichtian. We perform this temporal averaging because the paleogeographies of the Cenomanian and Maastrichtian are not reconstructed with a temporal resolution higher than a few million years. It is thus not possible to attribute a precise age to our Cenomanian (or Maastrichtian) paleogeography, which could equally appropriately represent a 97 or a 92 Ma paleogeography.

The Cenomanian is characterized by Atlantic and southern Indian Ocean ε_Nd values that range mainly between $\sim -\mathrm{5}$ and $\sim -\mathrm{6}$ in the intermediate ocean and $\sim -\mathrm{6}$ and $\sim -\mathrm{8}$ in the deep (Fig. 13). Exceptions to this are the anomalously low ε_Nd values recorded in the intermediate western equatorial Atlantic (Demerara Rise, MacLeod et al., 2008, 2011; Jiménez Berrocoso et al., 2010; Martin et al., 2012). The tropical Pacific has a high ε_Nd signature of $\sim -\mathrm{3}$; however, it is only represented by a single data point at Shatsky Rise (Murphy and Thomas, 2012).

Figure 13 Cenomanian and Maastrichtian ε_Nd compilation modified from Moiroud et al. (2016), with few additions (Tables S3 and S4). The ε_Nd values at each site are averaged between 100 and 90 Ma for the Cenomanian and between 75 and 65 Ma for the Maastrichtian. Site numbers are shown for clarity.

From the Cenomanian to the Maastrichtian, ε_Nd values generally decrease by ∼2 to 3 in the Atlantic and southern Indian oceans. In the Pacific Ocean, Maastrichtian ε_Nd values are $\sim -\mathrm{3.5}$ to −5.5 (Fig. 13). These ε_Nd trends have been the focus of numerous hypotheses suggesting the reorganization of ocean circulation through the Late Cretaceous (e.g., Robinson et al., 2010; MacLeod et al., 2011; Martin et al., 2012; Robinson and Vance, 2012; Murphy and Thomas, 2013; Voigt et al., 2013; Donnadieu et al., 2016; Moiroud et al., 2016). It has been suggested that the subsidence of large volcanic provinces, such as Kerguelen Plateau, could have decreased the supply of radiogenic material to the Southern Ocean and could have shifted the signature of Maastrichtian deep water masses formed in the South Atlantic (Robinson et al., 2010) or southern Indian Ocean (Murphy and Thomas, 2012) to lower values. Similar shifts toward lower ε_Nd values in the North Atlantic support hypotheses that suggest that by the Maastrichtian, the Central and South Atlantic had deepened enough to allow northward export of deep waters from the Southern Ocean to the North Atlantic (Robinson et al., 2010; Murphy and Thomas, 2012; Robinson and Vance, 2012).

The cessation of Pacific deep-water supply across the Caribbean Seaway in combination with an increased deep-water formation in the Atlantic and Indian sectors of the Southern Ocean has also been proposed to the ε_Nd shifts toward lower values (MacLeod et al., 2008; Donnadieu et al., 2016). Alternatively, these shifts could be explained by initiation of deep-water formation in the North Atlantic and invasion of the Southern Ocean by North Atlantic deep waters flowing across the equatorial Atlantic (MacLeod et al., 2005, 2011; Martin et al., 2012).

All of these hypotheses explain the similarity in deep-water ε_Nd values between the North Atlantic, South Atlantic, and Indian oceans in the Maastrichtian (Fig. 13) by greater communication between the basins (Robinson and Vance, 2012; Murphy and Thomas, 2013; Moiroud et al., 2016). Other records instead suggest that bathymetric barriers of the Rio Grande Rise (RGR)–Walvis Ridge (WR) system in the South Atlantic prevented deep north–south flow between the North Atlantic and the Southern Ocean until the Paleogene (Voigt et al., 2013; Batenburg et al., 2018), although recent work suggests that deep channels existed through the RGR–WR system in the Late Cretaceous (Moiroud et al., 2016; Pérez-Díaz and Eagles, 2017).

The opening of the Atlantic and Southern Ocean nonetheless played a major role in the convergent evolution of ε_Nd values in the Late Cretaceous by affecting intermediate and deep flow patterns as well as the residence time of water masses and, hence, local ε_Nd inputs such as boundary exchange.

5.2.2 Cenomanian circulation

In contrast to the model simulations of Donnadieu et al. (2016) and the observational hypotheses of Murphy and Thomas (2012, 2013) and Robinson et al. (2010), our Cenomanian simulation indicates deep-water formation in the southwestern Pacific, along the eastern coast of Australia, rather than in the South Atlantic or southern Indian Ocean (Fig. 3e). However, the deep-water pathway simulated in our Cenomanian simulation, with waters traveling from their southwestern Pacific source into the southern Indian and South Atlantic oceans following a strong westward current around the Australian continent, is reasonably consistent with existing ε_Nd proxy records. These deep waters would potentially have carried low ε_Nd values into the Indian and Atlantic sectors of the Southern Ocean because modern ε_Nd values of the margins close to the deep-water formation region in our Cenomanian simulation (eastern coast of Australia and Antarctic coast west of the Ross Sea) are typically between $\sim -\mathrm{7}$ and $\sim -\mathrm{20}$ (Jeandel et al., 2007; Roy et al., 2007). In the South Atlantic and southern Indian Ocean, deep-water ε_Nd values may have been modified by the addition of radiogenic contributions from recently active volcanic provinces (e.g., Kerguelen Plateau) that would raise the seawater value. Alternatively, it is possible that bathymetric barriers limited southwestern Pacific deep-water advection to the South Atlantic and southern Indian Ocean sufficiently to allow the ε_Nd signature of these deep water to be overprinted by regional ε_Nd supply in the Southern Ocean.

South Atlantic and southern Indian Ocean intermediate and deep sites do show a relatively large range of ε_Nd values (between $\sim -\mathrm{5}$ and $\sim -\mathrm{10}$, Fig. 13), and there is a wide range of possible ε_Nd sources with very different ε_Nd values. The African craton and Brazilian shield in the South Atlantic are unradiogenic (ε_Nd values < −10) (Jeandel et al., 2007), as are Antarctic terranes in the Atlantic and Indian sectors of the Southern Ocean (Roy et al., 2007). In contrast, the volcanic provinces of Walvis Ridge and Rio Grande Rise (O'Connor and Duncan, 1990; Murphy and Thomas, 2013; Voigt et al., 2013) and large igneous provinces of the Kerguelen Plateau and Rajmahal Traps (Mahoney et al., 1995; Coffin et al., 2002) exhibit more radiogenic values (ε_Nd values > −5). Precisely attributing the contribution of each source, including input of southwestern Pacific deep waters, to the South Atlantic and southern Indian Ocean ε_Nd values is, therefore, difficult.

Our Cenomanian simulation predicts an inflow of intermediate and deep waters into the North Atlantic from the Tethys and Mediterranean sections (Table S1 and Figs. 5a and 6a). These intermediate and deep waters mostly originate from the equatorial and tropical Pacific via an intense eastward current existing between ∼800 and 2400 m at the southern tip of Asia, which subsequently follows the eastern coast of Africa into the Neotethyan Ocean and the North Atlantic (Figs. 14a and 7c). Records from the equatorial Pacific (Murphy and Thomas, 2012) shows moderately high ε_Nd values (> −6) from the Cenomanian onwards. In addition, modern compilations of the ε_Nd signature of the continental margins in the eastern Mediterranean Sea (Ayache et al., 2016) and on the northeastern coast of Africa (Jeandel et al., 2007) indicate relatively radiogenic ε_Nd values (> −6). Inputs of radiogenic intermediate and deep waters from the Pacific into the North Atlantic via this Neotethyan pathway, regardless of whether sediment/water exchange in the Neotethyan Ocean may have contributed to their isotopic composition, provides a possible explanation for the ε_Nd signature of the deep North Atlantic (Fig. 13), which has more radiogenic values than the nearby North American and North African continents (Jeandel et al., 2007).

Figure 14 Cenomanian ocean circulation in (a) the northern Neotethyan Ocean at 1500 m depth, (b) the North Atlantic Ocean between 0 and 500 m depth, and (c) the North Atlantic Ocean between 1500 and 2000 m depth. Orange contours represent major pathways of water masses. Purple contours are the maximum winter MLD (500 m contours).

Intermediate and deep-water advection through the Neotethyan Ocean constitutes an alternative possibility to the direct deep-water advection from the Pacific to the North Atlantic through the Caribbean Seaway suggested by Donnadieu et al. (2016), which is problematic given that the Caribbean Seaway was probably closed to intermediate and deep-water flow as early as the Cenomanian (e.g., Buchs et al., 2018). However, other events may also have contributed to raising the ε_Nd values of North Atlantic intermediate and deep waters. In particular, volcanism related to the initial emplacement of the CLIP in the Caribbean Seaway during the Cenomanian could have supplied radiogenic material to the North Atlantic without requiring intermediate and deep-water exchange across the Caribbean Seaway or the Neotethyan Ocean. This input would raise the ε_Nd values of North Atlantic waters and could account for the high ε_Nd values ($\sim -\mathrm{5}$) observed in Cenomanian samples at Blake Nose in the intermediate North Atlantic (MacLeod et al., 2008). Another possible explanation for Blake Nose and other intermediate North Atlantic ε_Nd values could be a local supply of Pacific surface waters in the North Atlantic following a proto-Gulf Stream (Fig. 14b). The radiogenic surface signal could then have been transported to intermediate waters (Fig. 14c) via intermediate water formation in the North Atlantic (Figs. 3f and S5).

As pointed out in many studies, Demerara Rise and Cape Verde ε_Nd signatures stand out relative to other intermediate and deep sites (MacLeod et al., 2008, 2011; Jiménez Berrocoso et al., 2010; Martin et al., 2012). As in the simulation of Donnadieu et al. (2016), our Cenomanian simulation does not produce low-latitude intermediate or deep-water formation at Demerara Rise, as has been suggested by previous work (Friedrich et al., 2008; MacLeod et al., 2008, 2011; Martin et al., 2012). It does, however, show that Demerara Rise is bathed by a mixture of intermediate waters formed in the North Atlantic and originating from the Neotethyan Ocean, while the deeper Cape Verde site is mostly influenced by deeper waters from the Neotethys (Fig. 14c). It has been suggested that the low ε_Nd values at Demerara Rise could be due to boundary exchange with detrital material with extremely unradiogenic signature from the nearby Guyana shield (Donnadieu et al., 2016), possibly in conjunction with very restricted local circulation (Moiroud et al., 2016). Our model results support boundary exchange as an explanation for very low Demerara Rise values, but we cannot exclude the possibility that climate models are unable to reproduce low-latitude intermediate or deep-water formation at Demerara Rise because of missing processes or insufficiently detailed local paleogeography. Similarly, our results lead us to follow the suggestion that Cape Verde basin values could be driven by local boundary exchange close to the western African craton (Moiroud et al., 2016). We note that this conclusion is consistent with the results of Tachikawa et al. (1999, 2003), which report more unradiogenic values closer to the African continent at a site located in the high-organic-flux Mauritanian upwelling region rather than at a site located farther from the coast, which suggests a significant influence of boundary exchange processes in this region (Tachikawa et al., 2003).

5.2.3 Late Cretaceous circulation changes

The opening of the Atlantic and Southern oceans in our Maastrichtian simulations leads to an increased exchange of intermediate and deep waters between ocean basins (Figs. 5b and 6b), in line with previous model simulations (Donnadieu et al., 2016) and proxy-based evidence (e.g., Robinson et al., 2010; MacLeod et al., 2011; Friedrich et al., 2012; Martin et al., 2012; Robinson and Vance, 2012; Murphy and Thomas, 2013; Huber et al., 2018).

The evolution of the ocean circulation between the Cenomanian and the baseline Maastrichtian, 2× CO₂ Maastrichtian, or Deep Labrador Seaway experiments is reasonably consistent with the ε_Nd evolution to lower values. Because the 2× CO₂ Maastrichtian and Deep Labrador Seaway circulations are nearly identical to that of the baseline Maastrichtian experiment, we focus on the baseline Maastrichtian simulation. This simulation estimates higher rates of deep-water export from the southwestern Pacific to the Indian and Atlantic sectors of the Southern Ocean than the Cenomanian simulation (Fig. 6a–b). The absence of major changes in the provenance of deep currents between our Cenomanian and Maastrichtian model runs in the southern Indian and South Atlantic oceans suggests that the main cause of the observed decrease in ε_Nd in these basins might have been higher inputs of unradiogenic deep waters into the southern Indian and South Atlantic oceans driven by higher deep-water export rates and, therefore, less time for reactions with more radiogenic sediments (e.g., Haynes et al., 2020). Alternatively, the observed ε_Nd trend might be caused by the progressive subsidence of large igneous provinces, such as the Kerguelen Plateau, which would reduce the supply of radiogenic volcanic material to the Southern Ocean (Murphy and Thomas, 2013). These two hypotheses are not mutually exclusive and are difficult to test. However, we note that the shift toward lower ε_Nd values in the Indian and South Atlantic oceans is predicted by the slightly enhanced intensity of ocean circulation in the 2× CO₂ Maastrichtian simulation relative to the baseline Maastrichtian and is consistent with observational and model-based evidence for lower atmospheric CO₂ during the Maastrichtian (e.g., Wang et al., 2014; Tabor et al., 2016; Foster et al., 2017).

In our baseline Maastrichtian simulation, northward-flowing deep waters from the Southern Ocean dominate the Atlantic and could, therefore, advect low ε_Nd values to the North Atlantic and explain the observed ε_Nd signature shift in this basin (Figs. 6b and 13). This idea is consistent with previous arguments for the onset of an input of southern water masses into the North Atlantic (Robinson et al., 2010; Robinson and Vance, 2012; Murphy and Thomas, 2013). Indeed, in contrast to the separating role of the RGR–WR system on deep water masses suggested by Voigt et al. (2013) and Batenburg et al. (2018), gaps in the RGR–WR system in our Maastrichtian simulation are deep enough to allow northward flow of deep water.

Other studies have suggested that intermediate and deep waters could be sourced from high- (MacLeod et al., 2011; Martin et al., 2012) or low- (Friedrich et al., 2008; MacLeod et al., 2008, 2011) latitude regions in the North Atlantic, but deep-water formation there is not supported in our Maastrichtian simulation or in other recent coupled climate model simulations of the Late Cretaceous (Donnadieu et al., 2016; Lunt et al., 2016; Niezgodzki et al., 2017, 2019; Farnsworth et al., 2019). However, North Atlantic deep-water formation in the Cenozoic has been shown to be sensitive to details of North Atlantic configuration and bathymetry (Stärz et al., 2017; Vahlenkamp et al., 2018; Hutchinson et al., 2019). It is, therefore, possible that existing Late Cretaceous paleogeographic reconstructions are not sufficiently detailed, thereby inhibiting the modeled onset of North Atlantic deep-water production.

The Deep Caribbean Seaway and Deep Drake Passage simulations produce Pacific intermediate and deep waters that invade the Atlantic Ocean via northern or southern routes, respectively (Figs. 5d–e and 6d–e). This increased supply of Pacific waters into the Atlantic would be expected to increase the ε_Nd signature of the Atlantic basin, which is at odds with the observed ε_Nd decrease by ∼2 to 3 units from the Cenomanian to the Maastrichtian (e.g., Robinson et al., 2010; MacLeod et al., 2011; Martin et al., 2012; Robinson and Vance, 2012; Murphy and Thomas, 2013; Moiroud et al., 2016). Our simulations, therefore, argue against the presence of these deep gateways during the latest Cretaceous, in agreement with recent progress in the understanding of the geological history of these gateways but in notable contrast to the simulations of Donnadieu et al. (2016).

In the Deep Neotethys simulation as in the baseline Maastrichtian simulation, high volumetric flow rates of deep waters are exported from the southwestern Pacific to the Indian sector of the Southern Ocean (Fig. 6f), which, in conjunction with the subsidence of volcanic provinces, could explain the ε_Nd decrease in this basin. Because the Neotethyan Ocean is open to intermediate and deep circulation in this experiment, the deep North Atlantic is filled with westward-flowing deep waters from the Neotethyan Ocean, which then flow southward into the South Atlantic. These deep waters are composed of a mixture of southwestern Pacific deep waters with low ε_Nd values traveling across the Indian Ocean and of deep waters that have circulated in the tropical and equatorial Pacific Ocean and had their ε_Nd signature shifted toward higher values (e.g., Hague et al., 2012; Thomas et al., 2014; Haynes et al., 2020), before flowing into the eastern Neotethyan Ocean following the southern tip of Asia between ∼2000 and 3000 m (Figs. 15 and S14, Deep Neotethys Indonesian section). The low ε_Nd values observed in the Maastrichtian Atlantic could be consistent with a Deep Neotethys Seaway scenario if deep waters flowing into the North Atlantic were composed of a greater proportion of Pacific deep waters that traveled along the Indian Ocean and retained lower ε_Nd values than Pacific deep waters that traveled along the southern tip of Asia and acquired higher ε_Nd values. However, this hypothesis is less elegant and conceptually more complicated than the invasion of the North Atlantic by deep waters from the Southern Ocean with low ε_Nd values into the North Atlantic, as suggested by our baseline Maastrichtian (and Deep Labrador Seaway and 2× CO₂ Maastrichtian) simulation. In addition, the Deep Neotethys hypothesis is not easily reconciled with the geological context of a progressively resorbing Neotethyan Ocean during the Late Cretaceous (Stampfli, 2000).

Our Maastrichtian simulations offer no better solution to the low ε_Nd signature of Demerara Rise and Cape Verde records (MacLeod et al., 2008, 2011; Jiménez Berrocoso et al., 2010; Martin et al., 2012) than local boundary exchange processes within restricted basins (Donnadieu et al., 2016; Moiroud et al., 2016; Batenburg et al., 2018), at least until the extreme end of the Maastrichtian when a convergence of ε_Nd values from Demerara Rise and other North Atlantic sites is observed (MacLeod et al., 2011). Likewise, our simulations do not provide a particular solution to the high ε_Nd values recorded in Newfoundland basin in the Maastrichtian North Atlantic (Fig. 13). Thus, we concur with the suggestion that local processes involving more radiogenic material might contribute to this signal (Robinson and Vance, 2012), possibly as early as the Cenomanian (Fig. 13).

Figure 15 Deep Neotethys Seaway deep ocean circulation in the northern Indian and Neotethys oceans at (a) 3000 m depth and (b) 2500 m depth. Orange contours represent major pathways of water masses.

Figure 16 Surface, intermediate and deep ocean temperature difference (^∘C) between the 2× CO₂ Maastrichtian and Cenomanian simulations.

5.2.4 Oxygen and carbon isotopes

In contrast to the limited impact of specific gateway configurations on ocean temperatures in the Maastrichtian (Fig. 12a–d), regional surface and deep temperature changes of as much as ∼3 ^∘C are simulated between the Cenomanian and the Maastrichtian (Fig. 9a and S6). However, even in basins where these regional temperature changes are consistent with the direction of δ¹⁸O changes between the Cenomanian and the Maastrichtian, they fall short of explaining the amplitude of δ¹⁸O change observed in the proxy records (Huber et al., 2018). For example, the ∼1 ‰ to 1.5 ‰ positive benthic δ¹⁸O trend observed at Blake Nose (Huber et al., 2002, 2018) could in part be explained by the ∼2 ^∘C cooling predicted by our model in the North Atlantic (Figs. 10 and S6) between the Cenomanian and the Maastrichtian, but the parallel positive planktic δ¹⁸O trend is not reproduced in our simulations (Fig. 9a). Similarly, in contrast to proxy observations, our model does not predict any significant temperature change at Exmouth Plateau in the southern Indian Ocean or in the deep equatorial Pacific (Ando et al., 2013; Falzoni et al., 2016). In the Atlantic sector of the Southern Ocean, our model predicts a small cooling, which is consistent with the δ¹⁸O proxy record in terms of direction of change but not in amplitude (Huber et al., 2018).

As demonstrated in detail by Tabor et al. (2016), accounting for lower Maastrichtian atmospheric CO₂ levels allows better consistency between model results of the Cenomanian and Maastrichtian and observations (Fig. 16). North Atlantic Blake Nose temperature decreases by more than 4 ^∘C at ∼1000 m depth if Maastrichtian CO₂ levels are reduced by a factor of 2 relative to Cenomanian levels, in agreement with the ∼1 ‰ to 1.5 ‰ benthic δ¹⁸O trend. This cooling is paralleled by a ∼2 to 3 ^∘C surface cooling in agreement with the planktic δ¹⁸O record. With a halving of CO₂, the model also predicts a cooling of ∼2 to 2.5 ^∘C both at the surface and intermediate depth in the Indian Ocean at Exmouth Plateau and in the deep equatorial Pacific at Shatsky Rise, as well as a more pronounced cooling > 4 ^∘C at the surface and ∼3 ^∘C in the intermediate and deep ocean in the Atlantic sector of the Southern Ocean (Fig. 16 and Tables S3 and S4). These simulated temperature changes are in better agreement with proxy records (Ando et al., 2013; Falzoni et al., 2016; Huber et al., 2018) than in the absence of CO₂-induced cooling, in particular for benthic records. The amplitude of change in planktic δ¹⁸O between the Cenomanian and Maastrichtian is indeed generally larger than that of the benthic δ¹⁸O (Huber et al., 2018).

Part of the mismatch between simulated temperature changes and δ¹⁸O records may also pertain to the fact that foraminiferal δ¹⁸O is a proxy for temperature and seawater δ¹⁸O. Foraminiferal δ¹⁸O values are generally converted to temperatures using the consensus value of −1 ‰ for mean ice-free seawater δ¹⁸O (Shackleton and Kennett, 1975; Pearson, 2012), but regional deviations from the global mean seawater δ¹⁸O can exert a strong control on the conversion of foraminiferal δ¹⁸O values to ocean temperatures, in particular in the upper ocean. The mid-Cretaceous simulations of Zhou et al. (2008) with the GENESIS-MOM coupled model indicate significant surface variability in seawater δ¹⁸O in spite of the absence of a river routing scheme. Because precipitation and runoff are depleted in δ¹⁸O relative to seawater, the upper ocean could exhibit lower seawater δ¹⁸O in regions of high precipitation and/or high runoff input, with a substantial impact on reconstructed ocean temperatures (Huber et al., 2018).

Alternatively, if significant polar ice sheets developed during the Late Cretaceous, which is unlikely during the Cenomanian based on recent observational and model studies (e.g., MacLeod et al., 2013; Ladant and Donnadieu, 2016) but is more debated for the cooler climates of the Maastrichtian (e.g., Miller et al., 1999; Bowman et al., 2013; Ladant and Donnadieu, 2016; Huber et al., 2018), mean seawater δ¹⁸O may have shifted toward higher values. A positive shift in seawater δ¹⁸O would have reduced the magnitude of seawater cooling required to explain the increasing values in foraminiferal δ¹⁸O through the Maastrichtian. However, latest reviews suggest that in the absence of direct evidence for ice sheet and synchronicity between indirect evidence, Cretaceous ice sheets might only have existed, if ever, as small ice sheets with limited impact on seawater δ¹⁸O (Huber et al., 2018).

Finally, we note that CO₂-induced cooling may play a role in explaining the Cenomanian to Maastrichtian decrease in vertical δ¹³C gradients (Huber et al., 2018) because the temperature dependence of metabolic rates in ocean planktonic communities may have increased surface to deep δ¹³C gradient in warmer climates (John et al., 2013), by promoting increased rates of primary productivity, thereby enhancing surface δ¹³C values, and/or increased remineralization of organic matter, which would enhance the ¹³C depletion in the ocean interior.

In summary, the comparison of model results to planktic and benthic δ¹⁸O records confirms that prescribing lower atmospheric CO₂ levels in the Maastrichtian configuration is necessary to reproduce the cooling trend observed in the data. However, the absence of changes in ocean circulation with decreasing CO₂ levels and the limited changes in temperature produced by the deepening of gateways compared to that produced by lower CO₂ levels indicate that changes in both atmospheric CO₂ and paleogeography, likely with a strong influence from the nature of ocean gateways, are needed to reconcile model results and different proxy data into an internally consistent picture of evolving ocean circulation across the Late Cretaceous.

6 Conclusions

Our CCSM4 earth system model simulations of the Cenomanian and Maastrichtian demonstrate significant reorganizations of the deep and intermediate ocean circulation during the Late Cretaceous, which are predominantly controlled by the configuration of major oceanic gateways. Our model predicts continuous deep-water formation in the southwestern Pacific in the Late Cretaceous but shows that the Cenomanian to Maastrichtian interval witnessed the transition from an essentially zonal ocean circulation during the Cenomanian to one with increased meridional water exchanges during the Maastrichtian. We show that the simulated ocean circulation compares reasonably well to global ε_Nd records and that the Caribbean Seaway and Drake Passage were likely restricted to shallow circulation in the Maastrichtian, in agreement with current paleobathymetric knowledge (e.g., Buchs et al., 2018). In contrast, our simulations cannot discriminate whether deep connections existed across the Neotethyan Ocean on the basis of the comparison with ε_Nd records.

We are more confident in interpreting large basin-scale ε_Nd trends (such as the Atlantic and Indian oceans ε_Nd decrease between the Cenomanian and Maastrichtian) than local ε_Nd values. However, our interpretation of the larger patterns in the ε_Nd records is limited by several factors. First, paleogeographic uncertainties require that we average ε_Nd values over long time intervals. We are, therefore, bound to miss higher-frequency climatic and oceanic variability, which might explain regional ε_Nd signatures. Second, most of the neodymium signatures are between $\sim -\mathrm{5}$ and $\sim -\mathrm{10}$, which are relatively “middle-of-the-road” values that could be explained by a large number of plausible, not mutually exclusive, scenarios. Third, spatial and temporal resolution of the data is low for important intervals; there is a real need for increased Cretaceous ε_Nd records in particular from the south(western) Pacific and from the Indian Ocean, regions which are critically under sampled. These issues notwithstanding, direct comparison between ε_Nd records and oceanic currents is a step forward to understanding the ocean circulation of the Late Cretaceous, with future advances likely requiring specific modeling of the water mass signature in ε_Nd (Arsouze et al., 2007; Sepulchre et al., 2014; Gu et al., 2019).

Ultimately, our work highlights the critical impact of gateway configurations in the Late Cretaceous oceanic evolution. The geologic history of major ocean gateways and the continuous deep-water formation in the South Pacific in our simulations suggest that the Late Cretaceous trend in ε_Nd values in the Atlantic and southern Indian oceans was caused by subsidence of volcanic provinces and opening of the Atlantic and Southern oceans rather than changes in deep-water formation areas and/or reversal of deep-water fluxes. However, other plausible scenarios consistent with Late Cretaceous ε_Nd values remain, and new studies combining proxy records, detailed paleogeographic reconstructions, and ε_Nd modeling will, therefore, be key to improving our understanding of Late Cretaceous climates.