Middle Miocene climate–carbon cycle dynamics: Keys for understanding future trends on a warmer Earth?

The Miocene Epoch (23.03–5.33 Ma) stands out as a prolonged interval of relative global warmth, characterized by generally smaller polar ice sheets and a reduced latitudinal temperature gradient, thus potentially providing an analogue for future conditions on a warmer Earth. This interval was also marked by several episodes of major global climate change associated with fluctuations in ice volume, reorganizations of ocean circulation, and variations in transhemispheric heat and moisture transfer. The interval ca. 18–12.7 Ma, which encompasses the warmest climate phase of the Neogene and the transition to a much colder mode with extended polar ice cover, is especially intriguing because it offers insights into radically different mean states of climate variability on a warmer Earth. However, the driving mechanisms of these fundamental climate shifts, in particular, the relative roles of external radiative forcing versus internal feedback processes, remain to date highly controversial.

The Miocene climate optimum (ca. 16.9–14.7 Ma) represents the warmest interval of the Neogene period (e.g., Shackleton and Kennett, 1975; Flower and Kennett, 1994; Zachos et al., 2008) with mean annual temperatures estimated as ~4–8 °C above preindustrial values (e.g., Pound et al., 2012; Goldner et al., 2014). This extended period of warmth ended with stepwise expansion of the East Antarctic Ice Sheet and global cooling between ca. 14.7 and ca. 13.8 Ma (e.g., Kennett, 1977; Flower and Kennett, 1994; Holbourn et al., 2005; Shevenell et al., 2004, 2008) during the middle Miocene climate transition. The inception of a long-lasting positive carbon isotope excursion, known as the Monterey excursion (ca. 16.7–13.5 Ma), close to the onset of the Miocene climate optimum and the apparent covariance between δ¹⁸O and δ¹³C during the Miocene climate optimum originally gave rise to the hypothesis that enhanced burial of organic carbon and CO₂ drawdown ultimately drove middle Miocene Antarctic glaciation and global cooling (Vincent and Berger, 1985; Woodruff and Savin, 1991; Flower and Kennett, 1994). Low-frequency δ¹³C oscillations (~1‰ in amplitude) characterizing the Monterey excursion display a distinctive ~400 k.y. rhythm, suggesting that carbon reservoir changes were strongly modulated by eccentricity forcing (Woodruff and Savin, 1991; Holbourn et al., 2007). However, no consensus has yet been reached concerning the role of the terrestrial and marine carbon reservoirs and the involvement of nutrient/weathering feedbacks in driving middle Miocene changes in the carbon cycle and climate (e.g., Diester-Haass et al., 2009; Ma et al., 2011; Sosdian et al., 2020).

High-resolution marine records spanning the Miocene climate optimum and middle Miocene climate transition are still relatively scarce, largely due to the difficulty in recovering complete and well-preserved sedimentary successions over this interval of major climate change. Marine sequences are often affected by carbonate dissolution on account of the highly fluctuating carbonate compensation depth during that period (e.g., Pälike et al., 2012; Kochhann et al., 2016), or they exhibit strong diagenesis as a result of deep burial within the sediment column. Over the last decades, the International Ocean Discovery Program (IODP) and its predecessors the Integrated Ocean Drilling Program (also IODP) and the Ocean Drilling Program (ODP) have recovered some outstanding, undisturbed sedimentary successions from the Pacific, Indian, and Southern Oceans that provide new insights into the sequence and pacing of events across the Miocene climate optimum and middle Miocene climate transition (e.g., Shevenell et al., 2004, 2008; Holbourn et al., 2005, 2007, 2014, 2015; Miller et al., 2017). The primary objective of this work was to synthesize published and new high-resolution benthic foraminiferal isotope records from ODP/IODP drilled cores in the Pacific, Indian, and Southern Oceans and X-ray fluorescence (XRF) scanning data from a deep Pacific site to monitor interocean gradients and to reconstruct circulation changes for a better understanding of processes driving long- and short-term climate variability on a warmer Earth.

ODP Site 751 (57°43.56′S, 79°48.89′E, water depth: 1634 m), located in the Raggatt Basin on the Southern Kerguelen Plateau, lies ~900 km south of the present-day polar front (Fig. 1). The sedimentary succession recovered at this site consists of Lower to Upper Miocene “diatom nannofossil oozes” and had sedimentation rates of 1.5–2 cm/k.y. The succession is interrupted by a major hiatus, which spans the interval ca. 16–12.5 Ma (Schlich et al., 1989).

ODP Site 1146 (19°27.40′N, 116°16.37′E, water depth: 2092 m) is situated on the midcontinental slope of the northern South China Sea and is today located at the northwestern edge of the Western Pacific Warm Pool (Fig. 1). During the middle Miocene, the connection between the South China Sea and the western Pacific Ocean was more open than it is today, as the Bashi Strait only formed in the late Miocene (Hall, 1998; Wang et al., 2000). The continuous Miocene sequence of carbonate-rich hemipelagic sediments recovered at this site had an average sedimentation rate of ~2–3 cm/k.y. (Supplemental Fig. S1¹). Detailed site location, core recovery, and lithological descriptions were provided in Wang et al. (2000).

ODP Site 1171 (48°30′S, 149°06.69′E, water depth: 2150 m) was drilled on the southern tip of the South Tasman Rise (Fig. 1). In the middle Miocene, the site was located at a paleolatitude of 55°S with a backtracked paleodepth of 1600 m (Exon et al., 2001). The middle Miocene sequence recovered at this site consists predominantly of nannofossil oozes deposited with average sedimentation rates of ~1.4 cm/k.y. (Shevenell et al., 2008).

ODP Site 1236 (21°21.54′S, 81°26.17′W, water depth: 1323 m) and ODP Site 1237 (16°0.42′S, 76°22.69′W, water depth: 3212 m) are located on the Nazca Ridge off Peru (Fig. 1). During the middle Miocene, both sites were located ~10° westward of their present positions within the subtropical gyre away from the coastal upwelling regions of Peru (Mix et al., 2003). The sediment recovered at these sites is composed of unlithified, white to pale brown nannofossil oozes with an average of ~95% calcium carbonate (Mix et al., 2003). Average sedimentation rates over the middle Miocene were ~0.5–2 cm/k.y. (Supplemental Fig. S1).

IODP Site U1337 (3°50.01′N, 123°12.35′W, water depth: 4463 m) and IODP Site U1338 (2°30.47′N, 117°58.18′W, water depth: 4200 m) are located in the eastern equatorial Pacific Ocean (Fig. 1). Both sites were located close to the equator during the middle Miocene at slightly shallower paleodepths of ~4000–3500 m (Site U1337) and 3500–3000 m (Site U1338) and within or slightly northward of the narrow zone of equatorial upwelling (Pälike et al., 2010). The continuous carbonate-rich successions recovered at both sites exhibit average sedimentation rates of ~1.5–5 cm/k.y. (Supplemental Fig. S1).

IODP Site U1443 (5°23.01′N, 90°21.7′E, water depth: 2925 m) was drilled close to the location of ODP Site 758 on the Ninetyeast Ridge in the eastern equatorial Indian Ocean (Fig. 1). A continuous, carbonate-rich sedimentary succession spanning the Miocene was recovered at this site (Clemens et al., 2016). Average sedimentation rates over the middle Miocene were ~0.2–1.4 cm/k.y. (Supplemental Fig. S1).

In the modern ocean, the North Atlantic Deep Water (NADW), which is formed in the Norwegian-Greenland and Labrador Seas, and the Antarctic Bottom Water (AABW), which originates from the Weddell and Ross Seas (e.g., Worthington, 1970; Talley, 1999), play a crucial role in the global Meridional Overturning Circulation (MOC). The NADW, characterized by high δ¹³C, flows southward across the Atlantic Ocean toward the Weddell Sea, where it becomes entrained into the circum-Antarctic Current and contributes to the formation of AABW. Unlike NADW, AABW includes a smaller component of surface waters and thus has lower δ¹³C (~0.4‰; Kroopnick, 1985; Mantyla and Reid, 1983). AABW flows northward in the abyssal parts of the Pacific and Indian Oceans and is finally mixed into the central deep water at low latitudes in both ocean basins. Due to the long residence time of deep water in the North Pacific, remineralization of organic matter results in low oxygen concentrations and the lowest δ¹³C values (~–0.5‰) in the global ocean.

At the Pacific Ocean sites selected for this study, there is a clear relationship between δ¹³C and the age of deep-water masses (Fig. 2). Site 1146, located in the South China Sea at a water depth of 2092 m, is characterized by the lowest δ¹³C values of ~–0.5‰ and a dissolved oxygen content below 100 μmol/kg (Fig. 2). This site is bathed by poorly ventilated North Pacific Intermediate Water and Pacific deep water that enter the South China Sea through the Bashi Strait (Wyrtki, 1961). At eastern equatorial Pacific Sites U1337 and U1338 (>4000 m water depth), the dissolved oxygen content remains below 100 μmol/kg, and δ¹³C values reach ~0.5‰, reflecting the influence of AABW at this water depth. Sites 1236 and 1237, located in the subtropical southeastern Pacific Ocean, are characterized by δ¹³C values of ~0.25‰ and a dissolved oxygen content above 100 μmol/kg (Fig. 2), typical for Central Pacific Deep Water, which is a mixture of Northern and Southern Hemisphere–sourced deep-water masses.

The bottom water at eastern equatorial Indian Ocean Site U1443 (2925 m water depth), with a δ¹³C of ~0.25‰ and oxygen content of ~150 μmol/kg (Fig. 2), is composed of Indian Deep Water (IDW), which occupies the depth range from 3800 m to ~1500 m in the equatorial and northern Indian Ocean. IDW is characterized by high salinities reaching maxima of 34.8 practical salinity units (psu) in the southwestern part of the Indian Ocean and 34.75 psu in the southeastern part, where its upper limit rises to a depth of 500 m (Tomczak and Godfrey, 2003). Its temperature, salinity, and oxygen properties in the high-salinity core are virtually identical with those of NADW in the Atlantic sector of the Southern Ocean, indicating that IDW is mainly of NADW origin and is not formed in the Southern Ocean, in contrast to AABW, which fills the Indian Ocean below 3800 m (Tomczak and Godfrey, 2003). The northward flow of IDW, concentrated along the western boundary of the Ninetyeast Ridge, penetrates into the Northern Hemisphere, where it is modified by mixing with thermocline water from above and upwelling of AABW from below. The shallower Sites U1447 and U1448 (~1400 and ~1100 m water depths), located in the western Andaman Sea, have δ¹³C of ~–0.25‰ and a dissolved oxygen content below 100 μmol/kg, which is typical for IDW modified by vertical mixing (Fig. 2). In the Miocene, however, the influence of NADW on the temperature, salinity, and oxygen properties of IDW was substantially reduced, as the Southern Ocean was the main deep-water source for both the Indian and Pacific Oceans (e.g., Holbourn et al., 2013a).

The bottom water at Site 751 in the Indian Ocean sector of the Southern Ocean and at Site 1171 in the Pacific sector of the Southern Ocean is mainly composed of Circumpolar Deep Water (CDW), which is the main contributor to the Antarctic Circumpolar Current (ACC; Emery and Meincke, 1986). CDW is a mixture of deep and bottom waters originating from the subpolar Southern Ocean (Weddell Sea) with large contributions of Indian and Pacific Intermediate Water and NADW (Broecker et al., 1985). It is enriched in dissolved inorganic carbon and, thus, plays a fundamental role in defining the sensitivity of the global carbon cycle to the overturning circulation (Sigman et al., 2010). This is due to the fact that CDW can be either upwelled to the surface through wind-driven divergence or downwelled through mixing with dense waters from the continental shelf, thus contributing to the formation of AABW (MacGilchrist et al., 2019).

We integrated published benthic foraminiferal isotope records from eight sites in the Pacific and Southern Oceans with a new benthic δ¹³C record from Site U1443 in the Indian Ocean. The middle Miocene interval from Site U1443 was sampled for stable isotope analysis along the revised shipboard splice in 2 cm depth intervals. The processing of U1443 samples and the preparation of foraminifers for stable isotope analysis followed methods described in Lübbers et al. (2019). We measured δ¹³C in the epifaunal species Cibicidoides wuellerstorfi and Cibicidoides mundulus, except in eight samples within cores U1443A-15H and U1443B-15H, where these species were extremely rare or absent. In these eight samples, δ¹³C was measured on Planulina renzi, Cibicidoides havanensis, or Osangularia culter. Detailed methods and applied corrections are described in Supplemental Text 1 and Supplemental Table S1 (see ^{footnote 1}). Samples were analyzed at the Leibniz Laboratory for Radiometric Dating and Stable Isotope Research, Christian-Albrechts-University, Kiel, Germany, with a Finnigan MAT-253 coupled to a Carbo-Kiel type IV device. Results were calibrated against the standards NBS 19 and IAEA 603 and are reported in relation to the Vienna Peedee belemnite (VPDB) δ¹³C scale. We measured 96 duplicate samples in the middle Miocene interval of Site U1443, which indicated a standard error of ±0.1‰ for δ¹³C measurements. New U1443 stable isotope data are available at the Data Publisher for Earth and Environmental Science: https://doi.org/10.1594/PANGAEA.922280.

We also present X-ray fluorescence (XRF)–scanning elemental data from ODP Site 1237. Detailed methods for XRF scanning were described in Holbourn et al. (2013a), and XRF-scanning data are available at https://doi.org/10.1594/PANGAEA.922280.

Building on previously published chronologies (Holbourn et al., 2014, 2015; Kochhann et al., 2016), we developed a revised composite age model for Sites U1337 and U1338, based on correlation of the benthic foraminiferal δ¹⁸O series at each site to computed variations of Earth’s orbit (eccentricity, obliquity, and precession from the La04 orbital solution; Laskar et al., 2004; see Supplemental Figs. S2A and S2B). At Site U1338, only minor changes were implemented, and the revised chronology is almost identical to that published by Holbourn et al. (2014). At Site U1337, we implemented two corrections to the revised splice following Wilkens et al. (2013; also personal commun., 2017): (1) Core U1337C-20X is now aligned with Core U1337D-36X by tying U1337D-36X-6, 37 cm (322.37 m below seafloor [mbsf], 356.35 m composite depth [mcd] revised) to U1337C-20X-2, 49 cm (319.29 mbsf, 356.35 mcd revised). This correction results in an extension of the splice by 140 cm. (2) We additionally removed a duplicate interval of 50 cm between U1337C-20X-7, 4–6 cm and U1337A, 36X-3, 132 cm (363.7 m mcd revised). For constructing the composite chronology spanning 18–15.97 Ma at Site U1337 and 15.97–12.7 Ma at Site U1338, we used the δ¹⁸O minimum at ca. 15.97 Ma as a tie point (Supplemental Table S2; Supplemental Figs. S2A and S2B).

As tuning target, we constructed an eccentricity-tilt composite (ET) with no phase shift and with equal weight of eccentricity (E) and obliquity (O), and a value of 0.3 for the precession parameter P (ET+0.3P). After experimenting with several eccentricity-obliquity-precession mixes (e.g., Supplemental Figs. S2A and S2B), we selected this target curve, since it showed most similarity to our δ¹⁸O series. We did not use a negative sign for the precession parameter, since we assumed that Southern Hemisphere summer insolation drove global ice volume and climate change on a middle Miocene glaciated Earth with a dominant, dynamic Antarctic ice sheet. We did not adjust our tuning for possible phase lags between δ¹⁸O and insolation forcing, since the response time of a Miocene Antarctic ice sheet is unknown. We initially followed a minimal tuning approach to avoid imparting artificial changes in sedimentation rates (e.g., Muller and MacDonald, 2000) and then recognized that ET+0.3P maxima closely matched δ¹⁸O minima between primary tuning tie points and subsequently adjusted our age model for small offsets.

Here, we coded the ET+0.3P maxima and corresponding δ¹⁸O minima driven by 100 k.y. eccentricity with the first three digits of their ages in the La04 solution (Fig. 3; Supplemental Table S2). For instance, the ET+0.3 maxima at 12.724, 12.745, and 12.766 Ma in the La04 solution are labeled as 127a, 127b, and 127c (Supplemental Figs. S2A and S2B; Supplemental Table S2). For intervals between 14.7 and 14.0 Ma and between 16.8 and 16.1 Ma, where high-frequency variability was mainly driven by obliquity, we followed a comparable approach and labeled ET+0.3 maxima with the first four digits of their age. For instance, the ET+0.3P maximum at 14.347 Ma in the La04 solution is labeled 1434 (Supplemental Figs. S2A and S2B). A complete list of ET and ET+03P maxima, including labels, is provided in Supplemental Table S2. The coding of ET+0.3P maxima and δ¹⁸O minima was also integrated within the numbering scheme of 405 k.y. cycles from Wade and Pälike (2004) in Supplemental Figures S2A and S2B. Stable isotope data with revised ages are available at https://doi.org/10.1594/PANGAEA.922280.

To enable meaningful comparison of benthic isotope records from Sites 751, 1146, 1236, and 1237, we adjusted their independently tuned age models (Holbourn et al., 2005, 2007, 2013a; Miller et al., 2017) to the U1338-U1337 composite chronology. For this, we used the benthic δ¹³C records, which exhibited the most distinctive correlative structures (Fig. 4). This mainly resulted in only minor changes, which did not substantially alter the previously published age models. However, the revision of the oldest part of the Site 1146 chronology resulted in a substantial extension of the stratigraphic interval represented at this site, in contrast to the original age model, based on extrapolation of the sedimentation rates in the lowermost part of the sedimentary succession. In addition, a small hiatus of ~20–30 k.y. was identified in Site 1146 at the end of the major δ¹⁸O increase (ca. 13.85 Ma), which coincides with the first peak of the ~400 k.y. carbon isotope maximum (CM), known as CM6. Updated age models for Sites 751, 1146, 1236, and 1237 are provided in Supplemental Table S3 (see ^{footnote 1}).

The temporal resolution of the middle Miocene stable isotope record at Indian Ocean Site U1443 ranges between 1 and 15 k.y., with a mean value of ~5 k.y. For the youngest interval (13.1–12.7 Ma), we used the published isotope data and age model of Lübbers et al. (2019). Below this interval, the age model was constrained by 14 primary tie points between the U1443 δ¹³C record and the orbitally tuned combined U1338-U1337 splice. These tie points ranged between 157.36 and 136.95 mcd, with ages from 17.33 to 12.73 Ma (Supplemental Table S3). Between 17.33 and 13.91 Ma, mean sedimentation rate was 0.3 cm/k.y., rising to 1.4 cm/k.y. between 13.91 and 13.54 Ma, and then decreasing to 0.6 cm/k.y. between 13.55 and 13.14 Ma (Supplemental Fig. S1).

A widespread carbonate minimum, associated with a prominent hiatus (Lavender unconformity, ca. 17–16 Ma) and an increase in silica deposition, was initially detected at Deep Sea Drilling Program (DSDP) sites through the central equatorial Pacific Ocean (Mayer et al., 1985). Highly resolved records subsequently revealed that this massive carbonate dissolution event coincided with decreases in benthic foraminiferal δ¹⁸O starting at ca. 16.9 Ma and in δ¹³C between ca. 16.9 and 16.7 Ma (Figs. 3 and 4 e.g., Holbourn et al., 2015). The δ¹⁸O and δ¹³C shifts are also traceable in high-resolution records from ODP Site 1090 in the Southern Ocean (Billups et al., 2004), ODP Site 1237 in the southeastern Pacific Ocean (Holbourn et al., 2007), IODP Sites U1335 and U1337 in the equatorial Pacific Ocean (Holbourn et al., 2015; Kochhann et al., 2016), and IODP Site U1443 (Kochhann et al., 2021) in the Indian Ocean, thus supporting the global nature of these events.

Three distinct phases of climate variability can be identified during the Miocene climate optimum: (1) an initial warm phase from ca. 16.9–16.1 Ma, characterized by sustained low mean δ¹⁸O values of ~1.2‰ at Site U1338 (Fig. 3). The inception of warming at ca. 16.9 Ma was coupled to a negative shift in δ¹³C, lasting ~250 k.y., which preceded the onset of the first δ¹³C maximum (CM1) within the Monterey excursion at ca. 16.7 Ma (Figs. 3 and 4,Holbourn et al., 2015). (2) A sustained δ¹⁸O increase (~0.3‰ at Site U1338) from 16 to 15.7 Ma, contemporaneous with a marked δ¹³C increase (CM3, one of the most prominent δ¹³C maxima of the Monterey excursion), suggests a transient glacial expansion and CO₂ decrease (Mi2 event of Miller et al., 1991). (3) A prolonged interval of high-amplitude δ¹⁸O fluctuations (~1‰), corresponding to hyperthermal events primarily paced by short eccentricity and precession, followed between 15.6 and 14.7 Ma (Fig. 3; Supplemental Fig. S4 [see ^{footnote 1}]). During this interval, benthic δ¹³C displayed the low-frequency, ~400 k.y. variability characteristic of the Monterey excursion (Figs. 3 and 4). However, δ¹⁸O decreases (hyperthermal events) were also associated with sharp, transient drops in δ¹³C (Figs. 3 and 5) and intense carbonate dissolution episodes in the deep ocean (Holbourn et al., 2014; Kochhann et al., 2016), which were paced by 100 k.y. eccentricity (Fig. 3; Supplemental Fig. S4). This high-amplitude climate variability suggests that the Antarctic ice sheet remained highly dynamic and susceptible to insolation forcing during the Miocene climate optimum, in particular, between 15.6 and 14.7 Ma.

During the middle Miocene climate transition (14.7–13.8 Ma), stepwise cooling was associated with glacial expansion and establishment of a permanent East Antarctic Ice Sheet (e.g., Flower and Kennett, 1994; Shevenell et al., 2008; Holbourn et al., 2005). Changes in orbital cadence appear to have been instrumental in forcing the ocean-climate system over critical thresholds during this interval of fundamental climate change (Fig. 3). An initial increase in δ¹⁸O (~0.3‰) occurred at ca. 14.7 Ma, coincident with a marked decrease in amplitude variability and shortening of the dominant cyclicity from 100 to 41 k.y., lasting until ca. 14.1 Ma during an extended period of high-amplitude variability in obliquity, congruent with low variability in 100 k.y. eccentricity. The second, more important δ¹⁸O increase (~1‰) at ca. 13.8 Ma, associated with major expansion of the East Antarctic Ice Sheet and deep-water cooling (e.g., Flower and Kennett, 1994; Shevenell et al., 2008; Holbourn et al., 2005), also coincided with a change in eccentricity cadence (from 405 to 100 k.y.) and a decrease in the amplitude of obliquity variations (Holbourn et al., 2005). A strong precessional signal is embedded over the major δ¹⁸O decrease during the main glacial expansion from ca. 13.9 to 13.8 Ma (Holbourn et al., 2014). A further δ¹⁸O increase (~0.5‰) at ca. 13.1 Ma was associated with dampening of δ¹⁸O cycles (Fig. 3). The lower-amplitude variability in δ¹⁸O after the middle Miocene climate transition signals the establishment of more stable, extended ice sheets over Antarctica (Lewis et al., 2008; Holbourn et al., 2013b, 2018). The distinct, stepwise δ¹⁸O increases at ca. 14.7, ca. 13.8, and ca. 13.1 Ma, which mark the base of oxygen isotope zones Mi2a, Mi3, and Mi4 of Miller et al. (1991, 2017), track global cooling linked to progressive expansion of the East Antarctic Ice Sheet. Attendant sea-level falls associated with these major δ¹⁸O steps were estimated to be ~50 m for Mi2a, ~30 m for Mi3, and 20–30 m for Mi4 (Miller et al., 2020).

Comparison of δ¹³C profiles from eastern equatorial Pacific Site U1337 (water depth: 4463 m) and Southern Ocean Site 751 (water depth: 1634 m) revealed that the marked divergence between Pacific abyssal water masses (Site U1337) and Southern Ocean intermediate water masses that prevailed between ca. 18 and 16.9 Ma disappeared with the onset of the Miocene climate optimum and Monterey excursion (Fig. 6). During the Miocene climate optimum, the δ¹³C gradient between intermediate water masses (Southern Ocean Sites 751 and 1171 vs. eastern tropical Pacific Ocean Site 1237) and the δ¹³C gradient between abyssal water masses (eastern equatorial Pacific Site U1337 vs. Indian Ocean Site U1443) remained consistently reduced (Figs. 6 and 7; Supplemental Figs. S3A and S3B). This oceanwide decrease in stratification suggests a less vigorous meridional overturning circulation in response to a reduced meridional temperature gradient driving weaker Southern Hemisphere westerlies during the Miocene climate optimum. Furthermore, the low δ¹³C signature of Site 1146 deviates markedly from that of Pacific central water (Site U1338) throughout the Miocene climate optimum (Fig. 7), indicating the prevailing influence of a poorly ventilated intermediate water mass with a low preformed δ¹³C in the tropical northwestern Pacific Ocean. The presence of this poorly ventilated water mass, likely originating from the North Pacific Ocean, supports the scenario of a weaker meridional overturning cell within the Pacific Ocean during the Miocene climate optimum.

Model simulations also support a fundamentally different mode and strength of meridional overturning circulation during the warm and nearly ice-free Miocene climate optimum (Herold et al., 2012). These simulations suggested that the South Atlantic and Indian Ocean gyres were weaker and the interocean exchange via the ACC was reduced, while the South Pacific and Indian Ocean gyres were connected due to a more southerly position of Africa and a wider Indonesian gateway. Furthermore, the overturning of the upper water column was less intense and restricted to shallower depths due to weaker westerlies, resulting in reduced formation of Antarctic Intermediate Water and Upper CDW. By contrast, the northward flow of Lower CDW or AABW intensified, filling large parts of the global ocean to water depths as shallow as 1500 m (Herold at al., 2012). These modeling results are compatible with a previous compilation of benthic δ¹³C records that indicate the presence of a uniform deep-water mass in the world’s ocean during the Miocene climate optimum (Cramer et al., 2009) and with our comparison of δ¹³C profiles in the Pacific, Indian, and Southern Oceans (Figs. 6 and 7; Supplemental Figs. S3A and S3B).

Dissolution of deep-ocean carbonates in the eastern equatorial Pacific during the Miocene climate optimum was first recognized as the “Lavender” seismic reflector (Mayer et al., 1985). This prolonged interval of carbonate dissolution was associated with shoaling of the carbonate compensation depth, from a water depth of ~4.7 to 4.2 km (Lyle, 2003; Pälike et al., 2012). Widespread carbonate dissolution in the equatorial Pacific Ocean was initially attributed to the fractionation of silica and calcium carbonate between the Pacific and Atlantic Ocean Basins, possibly related to the formation of a proto–North Atlantic Deep Water (Mayer et al., 1985). However, more recent, high-resolution carbonate content records at equatorial Pacific Ocean Sites U1335, U1336, U1337, and U1338 and at southeastern Pacific Ocean Site 1237 revealed that Miocene climate optimum dissolution cycles were in phase with decreases in benthic δ¹⁸O and δ¹³C, and these changes were paced by short eccentricity (Holbourn et al., 2013a; Kochhann et al., 2016). Enhanced deep-ocean carbonate dissolution during warmer phases (Kochhann et al., 2016) suggests a tight coupling between climate and carbon cycle, possibly related to the release of light carbon (¹²C) from terrestrial mobile reservoirs, such as peat and permafrost, at eccentricity maxima (e.g., Zachos et al., 2010; DeConto et al., 2012). Thus, warmer phases of the Miocene climate optimum, mainly those between the thermal maximum at 15.6 Ma and the end of the Miocene climate optimum interval, are comparable to hyperthermal events that occurred during the Paleogene (e.g., Zachos et al., 2010).

The XRF-scanner–derived manganese (Mn) record from southeastern Pacific Ocean Site 1237 exhibited periodic enrichments closely matching δ¹⁸O decreases across the Miocene climate optimum (Fig. 8). Manganese enrichment during episodes of deep-water warming (hyperthermal events) is not only expressed in the Mn elemental counts and log(Mn/Ca), which could be explained by carbonate dissolution, but also in log(Mn/K) and log(Mn/Ti), which normalize Mn concentrations against typical terrigenous elements that are not affected by dissolution. Manganese enrichments during the Miocene climate optimum with enrichment factors exceeding the crustal molar Mn/Ti of 0.176 by 10–90 times were also reported by Chun and Delaney (2006), based on total digestion and inductively coupled plasma–mass spectrometry (ICP-MS) measurements of individual samples. Manganese enrichments relative to crustal abundances are usually formed by Mn oxides (MnO₂), indicating well-oxygenated conditions at the sediment surface (Calvert and Pedersen, 1993, 1996). However, authigenic Mn carbonates (MnCO₃) also precipitate within suboxic pore waters, promoting Mn precipitation and substantially increased Mn concentrations in environments with enhanced organic carbon burial and strongly oxygen-depleted pore waters (Boyle, 1983; Gingele and Kasten, 1994; Pedersen and Price, 1982). For instance, early diagenetic precipitates of Mn carbonates occur across the Paleocene-Eocene thermal maximum at intermediate- to deep-water sites (Chun et al., 2010; Pälike et al., 2014) and across the Cenomanian–Turonian oceanic anoxic event (OAE) 2 (Dickens and Owen, 1993) and the Aptian OAE 1a (Renard et al., 2005, 2007).

The mixing between oxygen-depleted water masses with high Mn solubility and well-oxygenated water masses with low Mn solubility has been suggested as a mechanism for Mn enrichment in carbonate environments associated with black shales during sea-level highstands. Manganese carbonates form by early diagenetic impregnation of carbonate substrate with anoxic waters saturated with MnCO₃ along steep redox gradients within the water column (Force and Cannon, 1988). Manganese enrichments in dysoxic pore waters have also been associated with enhanced early diagenetic release of methane from organic-rich sediments during hyperthermal events (Dickens, 2011; Pälike et al., 2014; Renard et al., 2005, 2007), and they appear to be a typical feature of warm, acidic, and dysoxic deep-water masses. The synchronous occurrence of eccentricity-paced transient warming events with deep-water ocean acidification and deoxygenation suggests a common triggering mechanism. Intervals of elevated atmospheric pCO₂ levels during and after periods of volcanic carbon dioxide release into the atmosphere would provide favorable boundary conditions for the widespread formation of such warm, acidic, and oxygen-depleted deep-water masses during the Miocene climate optimum.

Stepwise intensification of the δ¹³C gradient between intermediate- and deep-water masses in the Pacific Ocean following the Miocene climate optimum indicates a progressive increase in stratification and a more vigorous meridional overturning circulation following middle Miocene global cooling and glacial expansion. The δ¹³C gradient between deep-water masses at southeastern Pacific Ocean Site 1237 (3212 m water depth) and eastern equatorial Pacific Ocean Site U1338 (4200 m water depth) shows a sustained increase after ca. 14.7 Ma, following the end of the Miocene climate optimum (Fig. 7). Southern Ocean water masses (Site 1171, water depth: 2150 m) exhibit a similar diverging trend, though somewhat later after 14.2 Ma (Supplemental Fig. 3A). After ca. 13.6 Ma, intermediate waters at Site 1236 (water depth: 1323 m) increasingly diverged from deep-water masses at Site 1237 (water depth: 3212 m) and Site U1338 (water depth: 4200 m). Deep-water masses in the eastern equatorial Pacific Ocean (Site U1338) and Indian Ocean (Site U1443) also show a marked divergence after ca. 13.6 Ma (Supplemental Fig. 3B).

Steepening of the δ¹³C gradient between intermediate-water masses (Site 1236) and deep-water masses (Site 1237) in the southeastern Pacific Ocean was previously attributed to more vigorous entrainment of Pacific Ocean central water into the wind-driven ocean circulation, due to increased formation of intermediate and deep waters in the Southern Ocean during the middle Miocene climate transition (Holbourn et al., 2013a). The marked decrease in the δ¹³C gradient between northwestern Pacific Ocean Site 1146 (water depth: 2092 m) and eastern equatorial Pacific Site U1338 after 13.9 Ma further suggests that the northwestern Pacific Ocean became influenced by better-ventilated Pacific Ocean central water following global cooling (Fig. 7). A recent study, based on sedimentary coatings and fish teeth neodymium isotope data at Site U1438 (4700 m water depth) in the Philippine Sea, provided supportive evidence for the spreading of a deep, southern-sourced water mass to 27°N during the middle Miocene climate transition (Kender et al., 2018). In the North Atlantic, sea-surface temperature reconstructions at two strategic locations also indicated that a major shift in ocean heat transport occurred during the middle Miocene climate transition, probably due to intensification of the Atlantic Meridional Overturning Circulation and/or strengthening of the North Atlantic Current (Super et al., 2020). These findings suggest that a global reorganization of ocean circulation patterns occurred during the middle Miocene climate transition, which probably was instrumental in driving middle Miocene climate evolution.

Global cooling during the middle Miocene climate transition was also associated with a major improvement in deep-water ventilation and carbonate preservation (Holbourn et al., 2007; Kochhann et al., 2016). At deeper locations in the Pacific and Indian Oceans, epifaunal foraminifers, which are mainly suspension feeders, increased markedly after 13.9 Ma, while high-productivity foraminifers became more abundant at shallower sites (Smart et al., 2007; Singh et al., 2012; Holbourn et al., 2004, 2013a). A distinct 100 k.y. cyclicity is superimposed on this long-term trend, which appears to be associated with changes in circulation driven by eccentricity-paced variations of the expanded East Antarctic Ice Sheet (Holbourn et al., 2005, 2013a). The deep ocean remained poorly ventilated during warmer periods, whereas intermediate- and deep-water formation increased during colder intervals, resulting in improved ventilation and carbonate preservation in the deep ocean.

Global warming during the Miocene climate optimum has previously been linked to the intense phase of late early–middle Miocene flood basalt volcanism that led to the emplacement of the Columbia River Basalt Group along the Pacific Northwest of the United States (Hodell and Woodruff, 1994). It was proposed that intense volcanic degassing, resulting in atmospheric CO₂ loading, triggered rapid warming at the onset of the Miocene climate optimum. However, discrepancies between the basal ages of the lavas and the onset of the Miocene climate optimum have led to questions about volcanic forcing as the trigger of Miocene warmth (e.g., Mahood and Benson, 2017). The duration of the main phase of volcanic activity, initially estimated to have lasted ~1–2 m.y. (e.g., Barry et al., 2013), was recently reduced to a shorter interval (Mahood and Benson, 2017; Kasbohm and Schoene, 2018). While there is consensus for placing the waning phase of the main eruptive event at ca. 16 Ma, the timing for its onset remains controversial (Mahood and Benson, 2017; Kasbohm and Schoene, 2018; Moore et al., 2018; Cahoon et al., 2020). The Steen Basalts, long considered the earliest lavas of the Columbia River Basalt Group, were recently assigned basal ages of ca. 16.64 Ma (Mahood and Benson, 2017) and ca. 16.7 Ma (Kasbohm and Schoene, 2018), thus postdating the rapid warming at ca. 16.9 Ma at the onset of the Miocene climate optimum (Holbourn et al., 2015). However, a more recent study reported an extended stratigraphic range with substantially older ⁴⁰Ar/³⁹Ar basal ages of ca. 17.23 Ma for the Picture Gorge Formation (Cahoon et al., 2020), previously considered to postdate the Steens Basalts within the Columbia River Basalt Group (e.g., Reidel et al., 2013). In view of this new dating, the link between volcanism and climate warming at the onset of the Miocene climate optimum remains a possibility. However, a further difficulty resides in the volume of CO₂ released during eruption of the Columbia River Basalt Group, which has been estimated as too low to account for the amplitude of warming recorded during the Miocene climate optimum (e.g., Mahood and Benson, 2017). Despite these controversies, the sustained low benthic foraminiferal δ¹⁸O (global warming) from ca. 16.9–16.1 Ma through the main phase of volcanism, the marked δ¹⁸O enrichment (cooling) at ca. 16 Ma coincident with the waning of volcanic activity, the high amplitude of δ¹⁸O and δ¹³C variability (Fig. 3), and the elevated oceanic dissolved inorganic carbon (DIC) content (Sosdian et al., 2020) during the Miocene climate optimum all suggest that eruptive episodes substantially influenced late early to middle Miocene climate evolution, even if a direct, causal link between climate warmth and volcanism has not been established.

Today, neritic and pelagic environments sequester comparable amounts of carbonate, whereas the less extended neritic zone produced and accumulated approximately an order of magnitude less carbonate during glacial sea-level lowstands (Milliman and Droxler, 1996). If we assume that the global accumulation of deep-sea carbonate remained largely independent of sea-level changes and that shelf carbonate accumulation increased with expanding shelf areas, the deep ocean must have become depleted in alkalinity and carbonate during the warmer-than-present Miocene climate optimum, which was characterized by reduced ice volume and sea level ~50 m above today’s level (Miller et al., 2011, 2020). In addition to these changes in basin-shelf fractionation, the oceans must have absorbed a significant fraction of the excess atmospheric CO₂ released at the onset of the Miocene climate optimum and during hyperthermal events, thus becoming more acidic, with major implications for the marine carbonate cycle. As the lysocline and carbonate compensation depth shoaled, CaCO₃ production and burial would have decreased in both shelf and open-ocean environments, and carbonate dissolution would have increased, in line with the “coral reef hypothesis” (Berger, 1982a, 1982b; Opdyke and Walker, 1992). This hypothesis was originally developed to explain Pleistocene glacial-interglacial differences in atmospheric CO₂. Sea-level rises during Miocene hyperthermal events associated with a labile dynamic East Antarctic Ice Sheet would have likely exacerbated these trends, in a manner similar to Pleistocene glacial terminations (e.g., Opdyke and Wilkinson, 1988; Opdyke and Walker, 1992; Milliman and Droxler 1996). Furthermore, some of the organic carbon and methane hydrates accumulated in deep-sea sediments during cooler phases may have released isotopically light organic carbon as carbon dioxide or methane to the ocean-atmosphere system (e.g., Mawbey and Lear, 2013). Thus, the dynamics of carbonate accumulation on expanding shelf seas may have been a crucial driver of climate variability during the Miocene climate optimum.

The strong imprint of the 405-k.y.-long eccentricity cycle on low-latitude climate proxy records has previously been interpreted in terms of a direct response to low-latitude insolation changes as main drivers of monsoon variability (e.g., Wang et al., 2010). Based on energy-balance climate models, Short et al. (1991) attributed the larger eccentricity component in the temperature response of tropical regions (in addition to the precessional and obliquity components of direct insolation) to the effect of the twice-yearly passage of the sun across the equator. Box model results further supported the inference that carbon and nutrient fluxes to the deep ocean from continental weathering modulate the burial ratio of carbonate to organic carbon, thus imparting eccentricity-paced variability to δ¹³C in marine DIC (Ma et al., 2011). Enhanced continental weathering and river discharge linked to an intensified tropical hydrological cycle at insolation maxima during high eccentricity increase riverine DIC and nutrient and carbon fluxes to the ocean, resulting in δ¹³C depletion of marine bicarbonate in both the surface and deep ocean, and in an increase of the marine organic carbon pool. Based on reconstructions of surface-ocean carbonate chemistry, Sosdian et al. (2020) recently proposed that elevated atmospheric CO₂ from volcanic sources, warmer climate, and higher sea level led to increased surface-ocean DIC concentrations during the Monterey excursion, whereas organic carbon burial mainly occurred on an increased area of drowned continental shelves. These authors suggested that, in contrast to the entire Monterey excursion, the individual δ¹³C cycles (CM events), paced by long eccentricity, were driven by insolation-forced climatic cooling at high latitudes, which accelerated marine overturning circulation and upwelling of nutrient-rich intermediate waters in the tropics. This would have resulted in a positive carbon cycle feedback process of increased marine productivity and carbon burial, reduced surface-ocean DIC, and reduced atmospheric CO₂, consistent with the original Monterey hypothesis (Vincent and Berger, 1985).

Furthermore, the unusually high amplitude of δ¹³C variability at the 100 and 405 k.y. eccentricity bands during the Monterey excursion may have been related to changes in the size of the marine DIC pool (Paillard and Donnadieu, 2014). Weathering fluxes of nutrients and organic carbon are reduced at eccentricity minima, as marine overturning circulation intensifies and as carbonate and organic carbon burial rates increase in low-latitude upwelling regions, resulting in a global δ¹³C increase. A phase lag by ~50 ± 20 k.y. of the 405 k.y. δ¹³C cycle to eccentricity during the middle Miocene (Holbourn et al., 2007) further suggests a delay caused by the complex transfer of energy from the precessional band into the eccentricity band in combination with the smoothing effect of the long residence time of carbon in the ocean of ~100 k.y. (Wigley, 1976; Imbrie and Imbrie, 1980; Pälike et al., 2006). This phase lag is even larger in long-term carbon cycle reservoir modeling experiments for the Eocene–Miocene interval, which also showed that δ¹³C profiles vary as a function of cyclic burial of organic carbon in the ocean (Kocken et al., 2019).

Our integrated benthic foraminiferal isotope records also revealed that transitions to new background climate states during the Miocene climate optimum and middle Miocene climate transition were coupled to major δ¹³C increases/decreases, suggesting climate forcing past critical thresholds through carbon transfer between global reservoirs (Figs. 3 and 9). The most prominent climate turning points were (1) the negative δ¹³C excursion coincident with climate warming at the onset of the Miocene climate optimum (ca. 16.9 Ma), which lasted 250 k.y. and preceded CM1 at the onset of the Monterey excursion; (2) the global cooling synchronous with the δ¹³C maximum CM3 from 16.1 to 15.7 Ma, preceding the 100 k.y. hyperthermal mode; (3) the cooling step synchronous with the δ¹³C maximum CM5a, which coincided with the change from 100 to 41 k.y. variability at 14.7 Ma; and (4) the major cooling step from ca. 13.9 to 13.8 Ma, synchronous with the onset of the final δ¹³C maximum of the Monterey excursion (CM6), leading to a new mode of climate variability.

The positive δ¹³C excursions CM1, CM3, CM5a, and CM6 were all preceded by protracted episodes of global warming (Figs. 3 and 9). The onset of the positive δ¹³C shift of CM1 followed the intense warming episode at ca. 16.9 Ma, which marks the beginning of the Miocene climate optimum (Fig. 3). The major δ¹³C increases at the onset of CM3 (ca. 16.1 Ma), CM5a (ca. 14.7 Ma), and CM6 (ca. 13.8 Ma) also followed intervals of global warmth (Fig. 9). Global warming has been linked to intensification of the tropical hydrological cycle (enhanced monsoonal rainfall and tropical weathering; e.g., Ma et al., 2011), resulting in long-term carbon and nutrient transfer to the tropical oceans. We speculate that productivity increased when climate cooled during the onset of CM1, CM3, CM5a, and CM6, coincident with decreasing orbital eccentricity at the 405 k.y. and/or 100 k.y. band. Decreasing eccentricity enhances the meridional temperature contrast, promoting increased deep-water and intermediate-water formation in the Southern Ocean, a more dynamic overturning circulation, and an intensified biological pump.

Productivity increases at the onset of δ¹³C increases were likely linked to a chain of feedback mechanisms, including drawdown of atmospheric CO₂, global cooling, expansion of Antarctic ice sheets, intensification of westerlies, and increased silicate weathering on exposed shelves, which acted as an additional CO₂ sink (Wan et al., 2017). This is to some extent comparable to the amplifying role of atmospheric pCO₂ as a feedback mechanism, which boosted climate change during late Pleistocene glacial-interglacial cycles. However, the amplitude of Miocene δ¹³C change, which is in the range of ~1‰ for the initial δ¹³C increase (Fig. 3), is considerably higher than that in the late Pleistocene. Furthermore, cooler climate states (high benthic δ¹⁸O) were associated with intensified westerlies and trade winds in the Southern Hemisphere and with high δ¹³C on a Miocene unipolar glaciated Earth (Figs. 3 and 5; Supplemental Fig. S4). This phasing differs markedly from the typical out-of-phase variability of late Pleistocene glacial-interglacial cycles, when shifts in westerlies drove an intensified ACC during warm stages (Toggweiler, 2009). Climate cycles in the Miocene also were generally more symmetrical and exhibited faster rates of change than late Pleistocene cycles characterized by stepwise cooling followed by rapid warming, implying that they were driven by different carbon cycling feedback mechanisms.

We integrated high-resolution benthic foraminiferal δ¹³C and δ¹⁸O records from eight sites in the Pacific, Indian, and Southern Oceans to assess the evolution of interocean gradients and ocean circulation during the late early to middle Miocene in order to better understand the role of the marine carbon cycle and low-latitude processes in driving climate development on a warmer Earth. Alignment to a common orbitally tuned age model enabled discrimination of global trends (ice volume, meridional overturning circulation, and carbon reservoirs) and regional features (ventilation and oxygenation of water masses). Intensification of the δ¹³C gradient between intermediate and abyssal water masses in the Pacific and Southern Oceans prior to and after the Miocene climate optimum indicates enhanced stratification, associated with a more vigorous meridional overturning circulation during cooler periods. By contrast, ocean stratification decreased and the meridional overturning circulation substantially weakened during the Miocene climate optimum. X-ray fluorescence scanning data from a deep Pacific Ocean site additionally revealed that eccentricity-paced transient warming events during the Miocene climate optimum coincided with intense episodes of deep-water acidification and deoxygenation. These results suggest that changes in shelf-basin carbonate fractionation and weathering-induced carbon and nutrient fluxes from tropical monsoonal regions to the ocean, affecting the sequestration efficiency of the biological pump, played a crucial role in driving climate dynamics during the transition from an almost ice-free Earth during the Miocene climate optimum to an increasingly glaciated mode with permanent ice cover after the middle Miocene climate transition. This long-term perspective highlights the crucial role of internal feedback mechanisms in amplifying climate change and provides insights into the response of the climate–carbon cycle system to changes in mean-state background conditions on a unipolar glaciated, warmer Earth.

This research used samples and data provided by the International Ocean Discovery Program (IODP) and was funded by the Deutsche Forschungsgemeinschaft (grants HO233/1, KU649/27, and KU649/36-1). K.G.D. Kochhann and K.M. Matsuzaki gratefully acknowledge receipt of a Ph.D. scholarship (CNPq grant no. 244926/2013-1) from the National Council for Scientific and Technological Development, Brazil (Kochhann), and a fellowship (grant no. 19KK0089) from the Japan Society for the Promotion of Science (Matsuzaki). We thank Ken Miller and an anonymous reviewer for thoughtful, constructive comments.