Pliocene evolution of the tropical Atlantic thermocline depth

It has been hypothesized that global temperature trends are tightly linked to tropical thermocline depth, and that thermocline shoaling played a crucial role in the intensification of late Pliocene Northern Hemisphere glaciation. The Pliocene thermocline evolution in the Pacific Ocean is well documented and supports this hypothesis, but thermocline records from the tropical Atlantic Ocean are limited. We present new planktonic foraminiferal $\mathrm{Mg}/\mathrm{Ca}$, δ¹⁸O, and δ¹³C records from the late Pliocene interval at Ocean Drilling Program Site 959 in the Eastern Equatorial Atlantic (EEA), which we use to reconstruct ocean temperatures and relative changes in salinity and thermocline depth. Data were generated using surface-dwelling Globigerinoides ruber and subsurface-dwelling Neogloboquadrina dutertrei. Reduced gradients between the surface and subsurface records indicate deepening of the EEA thermocline at the end of the mid-Piacenzian Warm Period (mPWP; ∼ 3.3–3.0 Ma). We connect our late Pliocene records to previously published early Pliocene δ¹⁸O data from Site 959 and compare these to the Site 1000 in the Caribbean Sea. Over the course of the Pliocene, thermocline changes in the EEA and Caribbean Sea follow similar patterns, with prominent step-wise thermocline deepening between ∼ 5.5 and 4.0 Ma and gradual shoaling up to the mPWP, followed by minor deepening at the end of the mPWP. The tropical thermocline depth evolution of the tropical Atlantic differs from the Pacific, which is characterized by gradual basin-wide shoaling across the Pliocene. These results potentially challenge the hypothesized link between tropical thermocline depth and global climate. The mechanisms behind the periodically divergent Pacific and Atlantic thermocline movements remain speculative. We suggest that they are related to basin geometry and heterogenous temperature evolutions in regions from where thermocline waters are sourced. A positive feedback loop between source region temperature and tropical cyclone activity may have amplified tropical thermocline adjustments.

1 Introduction

The Pliocene (5.33–2.58 Ma) is the most recent geological epoch with substantially higher greenhouse gas concentrations (de la Vega et al., 2020; Martínez-Botí et al., 2015; Stap et al., 2016) and elevated global surface temperatures (Brierley et al., 2009; Dowsett et al., 2016; McClymont et al., 2020) compared to preindustrial times. This makes it an interesting interval to study for its potential analogies with our future climate (Burke et al., 2018). Typical features of the warm Pliocene ocean are the poleward expansion of the tropical warm pools and a reduction of the zonal sea surface temperature (SST) gradient in the Pacific (Brierley et al., 2009; Dekens et al., 2008; Fedorov et al., 2013). Another key feature of the Pliocene ocean is the deep tropical Pacific thermocline that gradually shoaled towards the end of the Pliocene (Dekens et al., 2007; Ford et al., 2012, 2015; Steph et al., 2010), while the meridional and zonal SST gradient steepened (Fedorov et al., 2015). It has been suggested that this shoaling reached a critical threshold around 3 Ma and played an important role in the onset of Northern Hemisphere glaciation (Fedorov et al., 2006; Philander and Fedorov, 2003; Steph et al., 2010). Moreover, because there is general coherence between Pliocene changes in meridional, zonal, and vertical (i.e., thermocline depth) temperature gradients in the Pacific, Fedorov et al. (2015) suggested that these temperature gradients are somehow mechanistically linked.

Interestingly, while the thermocline shoaled across the tropical Pacific, δ¹⁸O and $\mathrm{Mg}/\mathrm{Ca}$ records of planktic foraminifera show that the Caribbean thermocline deepened during the early Pliocene (Steph et al., 2010). This was interpreted by Steph et al. (2010) as a local phenomenon. However, early Pliocene δ¹⁸O records from the Eastern Equatorial Atlantic (EEA) also show a thermocline deepening (Norris, 1998). This suggests that the entire tropical Atlantic may have had a different early Pliocene thermocline evolution from the Pacific. If so, the link between tropical thermocline depth and Pliocene global temperature trends may need to be reconsidered. Late Pliocene thermocline records from the EEA are currently unavailable, yet necessary for a complete evaluation of Pliocene tropical thermocline depth evolution in a context of global climate changes.

We present the first records of late Pliocene thermocline depth variability in the EEA. We generated δ¹⁸O, δ¹³C, and $\mathrm{Mg}/\mathrm{Ca}$ records on surface-dwelling Globigerinoides ruber and subsurface-dwelling Neogloboquadrina dutertrei from Ocean Drilling Program (ODP) Site 959 sediments. We use the offsets between these records to infer changes in thermocline and nutricline depth and relative salinity changes. We evaluate our new records against previously published thermocline depth records from the Atlantic and Pacific to evaluate global patterns of thermocline depth changes across the Pliocene.

2 Methods

2.1 Site and modern oceanographic setting

We used sediments recovered at ODP Site 959 during Leg 159 in the Gulf of Guinea at ∼ 160 km offshore Ghana and Ivory Coast (3.62^∘ N, 2.73^∘ W; 2090 m depth; Mascle et al., 1996, Fig. 1). At present, this region is characterized by a shallow thermocline compared to the west Atlantic (Fig. 2). Surface ocean circulation at Site 959 is characterized by the eastward-flowing geostrophic Guinea Current that originates from the North Equatorial Counter Current and the Canary Current (Fig. 1). The thermal structure of the local water column varies seasonally in response to monsoon-induced changes in upwelling intensity (Fig. 1b). In West Africa, the intertropical convergence zone (ITCZ) shifts from 5–10^∘ N in February to 15–20^∘ N in August, which forces strong southwesterly monsoonal winds over the Gulf of Guinea in boreal summer (Gu and Adler, 2004). A major upwelling event occurs along the northern coast in summer, which is likely forced by a combination of local wind stress, wind-induced eastward-propagating Kelvin waves, and intensification of the Guinea Current (Djakouré et al., 2017; Verstraete, 1992), which together raise the thermocline. A shorter and weaker coastal upwelling event typically occurs in winter (Wiafe and Nyadjro, 2015). South of Site 959, upwelling occurs at the Equator in response to equatorial divergence and along the African west coast due to persistent Ekman pumping (Bakun and Nelson, 1991; Fig. 1).

Figure 1 General surface currents and upwelling areas in the Eastern Equatorial Atlantic. Blue shading: coastal upwelling in boreal summer; pink shading: equatorial upwelling; orange shading: permanent coastal upwelling. NEC: North Equatorial Current; NECC: North Equatorial Counter Current; CC: Canary Current; GC: Guinea Current; SEC: South Equatorial Current; EUC: Equatorial Undercurrent; SECC: South Equatorial Counter Current; BCC: Benguela Coastal Current. Green stippled lines mark the annual range of the Intertropical Convergence Zone (ITCZ). Figure after Norris (1998), Wagner (1998), and Wiafe and Nyadjro (2015). Map generated with Ocean Data View (Schlitzer, 2019).

Figure 2 (a) Upper ocean temperature (shaded contours; Locarnini et al., 2013) and salinity (black contours; Zweng et al., 2013) profiles across the black stippled line (see map) between Site 1000 and Site 959. The white stippled line shows the 20 ^∘C isotherm depth. Cross section generated with Ocean Data View (Schlitzer, 2020). (b) Seasonal vertical temperature and salinity profiles in the Gulf of Guinea (2.5^∘ N, 3.5^∘ W; Locarnini et al., 2013; Zweng et al., 2013). Symbols (left panel) mark the 20 ^∘C isotherm depth, which varies seasonally between ∼ 40–60 m.

2.2 Thermocline depth reconstructions

The thermocline is typically defined as the depth at which the vertical temperature change is at its steepest. This depth can be approximated with the 20 ^∘C isotherm in the modern ocean, which is typically done in modeling studies, but different isotherms may be more appropriate in warmer climates (e.g., von der Heydt and Dijkstra, 2011). A common approach for qualitative reconstructions of past thermocline depth variability is by paired tracing of surface ocean and shallow subsurface ocean conditions. We trace surface ocean conditions with Globigerinoides ruber and subsurface conditions with Neogloboquadrina dutertrei, which typically reside in the upper ocean and at thermocline depths respectively (Anand et al., 2003; Rebotim et al., 2017; Steph et al., 2009). Figure 3 conceptually shows how thermocline depth determines tracer offsets between surface and subsurface dwellers (Δsurface–subsurface). We use foraminiferal Δsurface–subsurface δ¹⁸O, $\mathrm{Mg}/\mathrm{Ca}$, and δ¹³C (hereafter Δδ¹⁸O, $\mathrm{\Delta }\mathrm{Mg}/\mathrm{Ca}$ and Δδ¹³C) to trace changes in thermocline and nutricline depth, and paired Δδ¹⁸O–$\mathrm{\Delta }\mathrm{Mg}/\mathrm{Ca}$ to trace relative salinity changes.

Figure 3 Conceptual manifestation of the thermocline depth as a $\mathrm{\Delta }\mathrm{Mg}/\mathrm{Ca}$ Δδ¹⁸O signal. Depth differentiation of the G. ruber and N. dutertrei habitats (stippled arrows) on a vertical temperature gradient (black line) sets a δ¹⁸O and $\mathrm{Mg}/\mathrm{Ca}$ offset in exported foraminiferal tests (solid arrows). The magnitude of this offset changes in response to thermocline shoaling and deepening.

2.3 Stable isotope analysis

Samples were taken at 5–20 cm intervals from cores 959C-5H and 959C-6H between 35.77 and 48.73 revised meters composite depth (rMCD; splice by Vallé et al., 2016). The studied interval of nannofossil or foraminifer ooze (Mascle et al., 1996) was dated using oxygen isotope stratigraphy (van der Weijst et al., 2020) and spans 3.5–2.8 Ma. The preservation of planktonic foraminifera in the Pliocene of 959 is generally very good; a large proportion of the tests have a “glassy” appearance, and G. ruber often still has spines around its aperture (Fig. S1 in the Supplement). These features suggest that diagenetic calcite overgrowth, dissolution, and recrystallization are minimal (Edgar et al., 2015).

For each measurement, 50–60 foraminiferal tests (G. ruber white s.s. and N. dutertrei) were picked from the 250–355 µm size fraction. Cleaning of G. ruber was performed at Utrecht University and that of N. dutertrei at the Royal Netherlands Institute for Sea Research (NIOZ) following slightly different protocols to remove clay and coccoliths. After gently crushing the tests between two glass plates, G. ruber was ultrasonically cleaned with ethanol and rinsed with distilled water. Crushed tests of N. dutertrei were ultrasonically cleaned three times with methanol, two times with hot (95 ^∘C) 1 % $\mathrm{NaOH}/{\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{2}}$ solution, and rinsed with Milli-Q in between each step. Stable isotope analyses were performed at Utrecht University on a Thermo Finnigan GasBench-II carbonate preparation device coupled to a Thermo Finnigan Delta-V mass spectrometer. The international IAEA-CO-1 standard was measured to calibrate to the Vienna Pee Dee Belemnite (VPDB) and the in-house NAXOS standard to correct for drift and track analytical precision, which was better than 0.1 ‰ for both carbon and oxygen isotopes.

2.4 $\mathrm{Mg}/\mathrm{Ca}$ analysis and calibration

Sediment samples were wet-sieved with tap water and dried overnight at a low temperature. Per sample, between 20 and 60 specimens of the surface-dwelling G. ruber (white, s.s.) and the thermocline-dwelling N. dutertrei were picked form the 250–355 µm size fraction and cleaned following the protocol of Barker et al. (2003). Tests were gently crushed between two glass plates and ultrasonically treated with methanol (3×), hot 1 % $\mathrm{NaOH}/{\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{2}}$ solution at 95 ^∘C (2×), 0.001 M HNO₃ (1×), and rinsed with Milli-Q in between each step. Trace element analyses were performed at the Royal Netherlands Institute for Sea Research (NIOZ) on an Element-2 ICP-MS. Note that the four M2 samples (identified by peak benthic δ¹⁸O; van der Weijst et al., 2020) were also used in a separate analysis for which a large quantity of carbonate was needed. For this reason, specimens of G. ruber from samples 959C/6H-2/102-104, 959C/6H-2/107-109, 959C/6H-2/112-114, and 959C/6H-2/117-119 were combined into a “M2 bulk sample”, and only 959C/6H-2/102-104 could also be analyzed separately. We calculate ocean temperatures using $\mathrm{Mg}/\mathrm{Ca}$ values (mmol/mol).

Ambient seawater temperature is typically the dominant control on the ratio of magnesium to calcium in foraminiferal shells (e.g., Hönisch et al., 2013; Tierney and Malevich, 2019). At present, there is no standard approach to calibrate $\mathrm{Mg}/\mathrm{Ca}$ to temperature, which complicates site-to-site comparisons (McClymont et al., 2020). However, there is general consensus that calibrating $\mathrm{Mg}/\mathrm{Ca}$ from fossil foraminifera to paleotemperature requires corrections for preferential removal of Mg over Ca during dissolution of calcite in the water column and sediment (Dekens et al., 2002; Regenberg et al., 2006) and the effect of the changing $\mathrm{Mg}/\mathrm{Ca}$ of seawater ($\mathrm{Mg}/{\mathrm{Ca}}_{\mathrm{sw}}$) on geological timescales in relation to rates of oceanic crust production (Hardie, 1996). Moreover, most calibrations assume a constant partitioning coefficient (D_Mg) at a given temperature, i.e., a linear relationship between $\mathrm{Mg}/{\mathrm{Ca}}_{\mathrm{sw}}$ and calcite $\mathrm{Mg}/\mathrm{Ca}$ ($\mathrm{Mg}/{\mathrm{Ca}}_{\mathrm{calcite}}$), but experiments suggest that D_Mg changes in response to $\mathrm{Mg}/{\mathrm{Ca}}_{\mathrm{sw}}$ (De Nooijer et al., 2017; Evans et al., 2016; Evans and Müller, 2012; Ries, 2004). To improve calibrations, species-specific D_Mg sensitivities and past $\mathrm{Mg}/{\mathrm{Ca}}_{\mathrm{sw}}$ need better constraints.

Here, we apply a constant Pliocene $\mathrm{Mg}/{\mathrm{Ca}}_{\mathrm{sw}}$ correction to the species-specific core top calibrations from Dekens et al. (2002), in which a constant partitioning coefficient is assumed, and a dissolution correction is embedded based on species, ocean basin, and core depth:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\mathrm{Mg}/\mathrm{Ca}\left(\mathrm{calcite}\right)=& \text{B}\cdot {\displaystyle \frac{\text{Past}\mathrm{Mg}/{\mathrm{Ca}}_{\mathrm{sw}}}{\text{Present}\mathrm{Mg}/{\mathrm{Ca}}_{\mathrm{sw}}}}\\ & \cdot e^{\mathrm{A}\left(T-C\cdot \text{core depth}\right)},\end{array}\end{array}$$

The slope (A) is 0.09 for G. ruber and 0.08 for N. dutertrei, the intersect (B) is 0.38 for G. ruber and 0.60 for N. dutertrei, and the dissolution parameter (C) is 0.61 for G. ruber and 2.8 for N. dutertrei (Dekens et al., 2002). The present-day $\mathrm{Mg}/{\mathrm{Ca}}_{\mathrm{sw}}$ is 5.2 mmol mol⁻¹, and we use the Pliocene estimate of 4.3 mmol mol⁻¹ (Evans et al., 2016). The core depth at Site 959 is 2.09 km. Note that G. ruber is calibrated to SST and N. dutertrei to 50 m depth subsurface temperature respectively (Dekens et al., 2002).

3 Results

Late Pliocene $\mathrm{Mg}/\mathrm{Ca}$ SST estimates from surface-dwelling G. ruber vary between 26.5 and 29.0 ^∘C and $\mathrm{Mg}/\mathrm{Ca}$ subsurface temperatures from N. dutertrei between 21.0 and 25.4 ^∘C (Fig. 4). Surface δ¹⁸O varies largely between −2.18 ‰ and −1.31 ‰ and subsurface δ¹⁸O between −1.19 ‰ and −0.14 ‰. Surface δ¹⁸O and $\mathrm{Mg}/\mathrm{Ca}$ barely show a long-term trend, although both record distinct cooling during the M2 glacial and peak temperatures during the KM3 interglacial. Subsurface δ¹⁸O and $\mathrm{Mg}/\mathrm{Ca}$ both show a warming trend following the M2 glacial, although it is more gradual in $\mathrm{Mg}/\mathrm{Ca}$-SubST than in δ¹⁸O-SubST. Surface and subsurface δ¹³C values vary between 0.91 ‰–2.13 ‰ and 0.56 ‰–1.45 ‰ respectively. Both show an increasing trend in the older part of the interval, and abruptly drop around KM3. In the younger part of the interval, surface δ¹³C remains stable, while subsurface δ¹³C increases.

Figure 4 (a) Surface (G. ruber; blue) and subsurface (N. dutertrei; orange) δ¹⁸O and $\mathrm{Mg}/\mathrm{Ca}$ temperature. (b) Surface (G. ruber; green) and subsurface (N. dutertrei; pink) δ¹³C records. Bold lines in (a) and (b) were smoothed with a LOESS moving regression model using a span of 0.5 in PAST software (Hammer et al., 2001). (c) Δsurface–subsurface offsets from (a) and (b) with Site 959 benthic δ¹⁸O record (van der Weijst et al., 2020) as grey silhouette for reference. Vertical bands indicate the M2 glacial and KM3 interglacial stage. Note that planktic δ¹⁸O and temperature axes are evenly scaled according to the δ¹⁸O-T relationship of Shackleton (1974).

The Δsurface–subsurface $\mathrm{Mg}/\mathrm{Ca}$, δ¹⁸O, and δ¹³C records (Fig. 4) follow a similar trend towards smaller offsets. Both $\mathrm{\Delta }\mathrm{Mg}/\mathrm{Ca}$ and Δδ¹⁸O signal a vertical temperature gradient decrease from ∼ 5.0 to 3.5 ^∘C. The Δδ¹⁸O record decreases after the M2 glacial, thereby leading the $\mathrm{\Delta }\mathrm{Mg}/\mathrm{Ca}$ and Δδ¹³C decrease around the KM3 interglacial by ∼ 150 kyr. Moreover, Δδ¹⁸O stabilizes in the younger part of the interval, while $\mathrm{\Delta }\mathrm{Mg}/\mathrm{Ca}$ and Δδ¹³C continue to decrease.

4 Discussion

4.1 Late Pliocene thermocline, nutricline, and salinity changes in the Eastern Equatorial Atlantic

The decreasing surface–subsurface $\mathrm{Mg}/\mathrm{Ca}$ and δ¹⁸O gradients between 3.3 and 2.8 Ma at Site 959 reflect thermocline deepening in the Eastern Equatorial Atlantic (Fig. 4). These decreasing vertical gradients are driven by a gradual warming trend in the subsurface ocean: G. ruber surface ocean $\mathrm{Mg}/\mathrm{Ca}$ and δ¹⁸O records remain stable across the interval, whereas the N. dutertrei subsurface records show an increase in $\mathrm{Mg}/\mathrm{Ca}$ and decrease in δ¹⁸O following the M2 glacial. Crude δ¹⁸O-based temperature estimates calculated with the linear Shackleton (1974) calibration have similar values and amplitude of variation as the $\mathrm{Mg}/\mathrm{Ca}$ estimates, but the decreasing Δδ¹⁸O trend leads decreasing $\mathrm{\Delta }\mathrm{Mg}/\mathrm{Ca}$ (Fig. 4).

Besides analytical and calibration errors, several confounding factors limit the use of Δδ¹⁸O and $\mathrm{\Delta }\mathrm{Mg}/\mathrm{Ca}$ as quantitative proxies for thermocline depth, including vertical migration of foraminifera in the water column (e.g., Steph et al., 2009), variability in the mass and/or preservation of the high-δ¹⁸O and low-$\mathrm{Mg}/\mathrm{Ca}$ secondary calcite crust that N. dutertrei precipitates at the end of its life cycle (Jonkers et al., 2012; Steinhardt et al., 2015), seawater carbon chemistry (Spero et al., 1997), and light intensity (Spero and Michael, 1987). Moreover, a discrepancy is observed in depth habitat ranges of planktonic species based on tows versus those reconstructed based on their δ¹⁸O values in core tops (Rebotim et al., 2017). Crucially, however, Δδ¹⁸O and $\mathrm{\Delta }\mathrm{Mg}/\mathrm{Ca}$ perform well in predicting larger thermocline patterns in the tropical Atlantic (Anand et al., 2003; Steph et al., 2009).

In addition to ambient temperature, calcite δ¹⁸O is determined by local seawater δ¹⁸O(δw), which is in turn tightly linked to salinity (Bigg and Rohling, 2000). Recently, Dämmer et al. (2020) advised against using coupled $\mathrm{Mg}/\mathrm{Ca}$–δ¹⁸O data for salinity reconstructions, because cumulative uncertainties from aforementioned confounding factors are too large to reliably reconstruct measured δw. We are therefore reluctant to calculate absolute salinities with our Pliocene data. However, these uncertainties are mostly relevant when applying this method quantitatively to discrete samples. To reduce these uncertainties, we use a semi-quantitative approach to assess regional long-term salinity trends by calculating the relative vertical salinity gradient (Δsalinity) from the discrepancy between $\mathrm{\Delta }\mathrm{Mg}/\mathrm{Ca}$ and Δδ¹⁸O (Fig. 5). We assume that $\mathrm{\Delta }\mathrm{Mg}/\mathrm{Ca}$ best approaches Δtemperature, despite the salinity effect on $\mathrm{Mg}/\mathrm{Ca}$ (Gray et al., 2018; Hönisch et al., 2013), and that the residual Δδ¹⁸O is explained by salinity changes. Note that the intersects of the $\mathrm{\Delta }\mathrm{Mg}/\mathrm{Ca}$ and Δδ¹⁸O temperature reconstructions in Fig. 5 may not correspond to a vertical salinity gradient of 0 due to uncertainties in both the δ¹⁸O paleotemperature and the $\mathrm{Mg}/\mathrm{Ca}$ calibrations. Rather, the intercepts indicate an unknown baseline Δsalinity to which the rest of the record can be compared.

Figure 5 (a) Smoothed $\mathrm{Mg}/\mathrm{Ca}$ and δ¹⁸O records from Site 959. (b) Relative magnitude of the salinity gradient between surface and subsurface, calculated by deconvolving Δsurface–subsurface δ¹⁸O and $\mathrm{Mg}/\mathrm{Ca}$, and assuming a modern tropical Atlantic δ¹⁸O-S slope (0.11 ‰ per salinity unit; solid line) and a doubled slope (dashed line). (c) LOESS smoothed (span of 0.5 in PAST software, Hammer et al., 2001) thermocline and nutricline tracers from Site 959 (as in Fig. 4c) and Site 1000 in the Caribbean Sea (stippled lines, data from Steph, 2005). Site 1000 Δδ¹³C was steeper than at Site 959 and was corrected with −0.4 ‰ to aid visual comparison.

The magnitude of variability in the Δsalinity record is strongly affected by the assumption of the δw–salinity (S) slope (Fig. 5). In the modern tropical Pacific, the δw–S slope varies with a factor >3 (0.09 ‰–0.32 ‰ per salinity unit) between across the equatorial Pacific, because the precipitation-forced isotopic “amount effect” (lower δw with higher tropical precipitation rates) is stronger in the western than eastern Pacific (Conroy et al., 2017). In the tropical Atlantic, however, slopes of 0.11 ‰ per salinity unit (Craig and Gordon, 1965) and 0.15 ‰ per salinity unit (LeGrande and Schmidt, 2006) have been reported for the tropical Atlantic region, but the zonal variability of the δw–S relationship is poorly documented. Using the slope provided by Craig and Gordon (1965; 0.11 ‰ per salinity unit), we calculate that Δsalinity was reduced by ∼ 2 salinity units during the mid-Piacenzian Warm Period (mPWP; 3.264–3.025 Ma) compared to the bordering intervals (Fig. 5). At present, the vertical salinity gradient is ∼ 1 in the Gulf of Guinea (Fig. 1b), suggesting that the Pliocene baseline salinity gradient was stronger than at present and/or that the δw–S slope was steeper. A doubling of the assumed tropical Atlantic δw–S slope (0.22 ‰ per salinity unit) is needed to solve the Δsalinity peak to the modern gradient of ∼ 1. It falls within range of modern tropical Pacific slopes and may be realistic for the interglacial Pliocene Atlantic, which was characterized by intensified monsoonal precipitation over West Africa (Haywood et al., 2020).

In the absence of major circulation changes, the Δδ¹³C record mainly reflects changes in the biological carbon pump and nutricline depth, which in the modern tropical ocean is tightly coupled to thermocline depth (Wilson and Adamec, 2002). We here consider coupled thermocline and nutricline deepening to be the most likely explanation for the observed synchronous changes in $\mathrm{\Delta }\mathrm{Mg}/\mathrm{Ca}$ and Δδ¹³C around KM3. The relative lead of the Δδ¹⁸O decrease may be explained by subsurface freshening (i.e., deepening of the halocline, Fig. 2b) in response to increased precipitation rates during the mPWP. This implies that $\mathrm{\Delta }\mathrm{Mg}/\mathrm{Ca}$ may outperform Δδ¹⁸O as a proxy for thermocline depth, but salinity changes would have to be strong in order to fully obscure the thermocline depth signal in the Δδ¹⁸O record. Despite temporarily increased precipitation during the mPWP, the net Δtemperature estimates from Δδ¹⁸O and $\mathrm{\Delta }\mathrm{Mg}/\mathrm{Ca}$ are equal between 3.5 and 2.8 Ma.

4.2 Dissimilar evolutions of the tropical Atlantic and Pacific thermocline

Our results show that the Eastern Equatorial Atlantic thermocline and nutricline deepened following the warmest mPWP interglacials, while SST and deep ocean records show global cooling (Herbert et al., 2016; Lisiecki and Raymo, 2005). Late Pliocene thermocline and nutricline changes are similar at Site 959 in the Eastern Equatorial Atlantic and Site 1000 in the western tropical Atlantic (Fig. 5), indicating that tropical thermocline deepening occurred across the entire basin. To further explore global patterns of Pliocene tropical thermocline movements, we connect our new Site 959 δ¹⁸O data to the early Pliocene records of Norris (1998) in Fig. 6. Across the Pliocene, thermocline depth evolution was similar in the Eastern Equatorial Atlantic and Caribbean Sea (Fig. 6). The late Pliocene ∼ 0.4 ‰ Δδ¹⁸O decrease at Site 959 and Site 1000 is relatively small compared to the fluctuations in the latest Miocene and early Pliocene (Fig. 6). Both sites show a net ∼ 1.5 ‰ Δδ¹⁸O decrease between 5.5 and 4.0 Ma, indicating significant thermocline deepening. This is followed by an interval of moderate thermocline shoaling, before deepening again in the late Pliocene. With the exception of the late Pliocene interval at Site 1241 (Fig. 6), tropical Pacific proxy records reflect overall Pliocene thermocline shoaling (Ford et al., 2012, 2015; LaRiviere et al., 2012; Steph et al., 2006b, 2010). Although changes in Δδ¹⁸O may not always be proportional to vertical movement of the thermocline (see Sect. 4.1), the early Pliocene fluctuations in the Atlantic records are appreciably larger than the long-term shifts towards lower Δδ¹⁸O values in the Pacific (Fig. 6). This suggests that the tropical Atlantic thermocline underwent major changes during the early Pliocene, which were not paralleled in magnitude by late Pliocene changes in the Atlantic or by any Pliocene thermocline change in the Pacific.

Figure 6 Compilation of tropical Δδ¹⁸O data. (a) Eastern Equatorial Atlantic (EEA) Site 959 (this study and Norris, 1998). (b) Western Tropical Atlantic (WTA) Site 1000 (Steph, 2005, and Steph et al., 2006a). (c) Eastern Equatorial Pacific (EEP) Site 1241 (Steph et al., 2006b). (d) Western Equatorial Pacific (WEP) Site 806 (LaRiviere et al., 2012). (e) LOESS smoothed records (span of 0.5 in PAST software; Hammer et al., 2001). Surface Trilobatus sacculifer records were corrected with +0.33 ‰ to G. ruber based on their offset in the modern ocean (Steph et al., 2009). Blue square marks late Pliocene interval covered in Figs. 4 and 5.

It was suggested by Fedorov et al. (2015) that global temperature trends are tightly linked to tropical thermocline depth and that thermocline shoaling played a crucial role in the onset of Northern Hemisphere glaciations around 3 Ma (Fedorov et al., 2006). Heat that is gained in the tropics must be balanced by heat loss at high latitudes, and because more heat can be gained by an ocean with a shallow thermocline, tropical thermocline shoaling should lead to surface ocean cooling at high latitudes (Boccaletti et al., 2004; Fedorov et al., 2006; Philander and Fedorov, 2003). While this theory is supported by the majority of the Pacific thermocline records (Ford et al., 2012, 2015; LaRiviere et al., 2012; Steph et al., 2006b, 2010), it appears to be inconsistent with basin-wide deepening of the tropical Atlantic thermocline (Fig. 6). The link between ocean SST gradients, tropical thermocline depth, and ocean heat transport may be different in the Atlantic and Pacific oceans, because the asymmetric geometry of the Atlantic basin leads to a different pattern of tropical thermocline ventilation (Harper, 2000). Also, Atlantic Meridional Overturning Circulation (AMOC) leads to northward ocean heat transport in both hemispheres, in contrast to the Pacific Ocean (Forget and Ferreira, 2019). However, because the tropical Pacific is larger than the tropical Atlantic, Pacific thermocline shoaling may have been sufficient to balance heat loss at high latitudes, even if partly counteracted by Atlantic thermocline deepening.

4.3 Controls on Pliocene thermocline changes

4.3.1 Central American Seaway closure

It has been hypothesized that Central American Seaway (CAS) closure played a major role in the global tropical thermocline depth (Steph et al., 2010; Zhang et al., 2012). In this scenario, salt transport to the North Atlantic increased as a consequence of reduced inflow of relatively fresh Pacific surface waters, which in model simulations promotes the production of NADW (Lunt et al., 2007; Sepulchre et al., 2014; Steph et al., 2010; Zhang et al., 2012). Steph et al. (2010) reasoned that an increased volume of NADW was associated with a greater volume of the “cold water sphere”, which raised the tropical thermocline everywhere except for the Caribbean region. There, reduced inflow of cold Pacific waters caused the thermocline to deepen locally. However, our study shows that tropical thermocline deepening occurred across the basin (Fig. 6). Even if CAS closure had promoted NADW production and/or AMOC strength, it is not clear how this, in turn, would have affected tropical thermocline depth, as both negative (Lopes dos Santos et al., 2010) and positive (Venancio et al., 2018) relationships between thermocline depth and AMOC strength have been inferred from proxy data. Moreover, it is unclear if CAS closure could have had an opposite effect on tropical Atlantic and Pacific thermocline depths, as this conflicts with modeling studies (Steph et al., 2010). A dominant role of CAS closure in shaping the tropical thermocline evolution in the Atlantic and Pacific basins is not consistently supported by data and models. Future work could be aimed at modeling the spatial extent of the region where the inflow of cold Pacific waters is reduced in response to CAS closure. Moreover, it is vital to constrain the chronology of CAS closure (Montes et al., 2015; O'Dea et al., 2016) and determine during which intervals this process might have had extensive and long-lasting oceanographic effects (Molnar, 2008; Sepulchre et al., 2014).

Figure 7 Regionally sorted midlatitude temperature records. (a) North Atlantic; (b) South Atlantic; (c) North Pacific; (d) South Atlantic. The records from Site 590 and 516 are $\mathrm{Mg}/\mathrm{Ca}$-based (T. sacculifer) from Karas et al. (2017); all other records are ${\mathrm{U}}_{\mathrm{37}}^{{\mathrm{k}}^{\prime }}$ temperature records from the Herbert et al. (2016) compilation.

4.3.2 Source region SSTs and tropical cyclones

Tropical thermocline waters are sourced from midlatitude surface waters (Harper, 2000). It has therefore been suggested that extratropical cooling contributed to the Pliocene shoaling of the tropical Pacific thermocline (Ford et al., 2012, 2015). During the early Pliocene, midlatitude SSTs were indeed generally higher than preindustrial (Brierley et al., 2009). However, Pliocene midlatitude temperature evolutions were regionally variable (Herbert et al., 2016; Karas et al., 2017), with a notable asymmetry between the Northern and Southern Hemisphere (Pontes et al., 2020). At present, Pacific thermocline waters are sourced from midlatitude surface waters in both the Northern and Southern Hemisphere, whereas in the Atlantic, the thermocline is predominantly sourced from the Southern Hemisphere as a consequence of the asymmetric basin geometry (Harper, 2000). Midlatitude SST records from potential Pacific tropical thermocline source regions predominantly register cooling between 5.5 and 2.8 Ma. In contrast, South Atlantic Site 1088 registers warming during the latest Miocene–early Pliocene, and during the mPWP (Fig. 7), in tandem with tropical Atlantic thermocline deepening. However, South Atlantic cooling between ∼ 4.5 and 3.5 Ma (Fig. 7) is not mirrored by thermocline shoaling, suggesting that source region SST was not likely the only driver of tropical Atlantic thermocline depth.

Tropical thermocline depth is also affected by cyclone activity, which, in turn, is linked to the latitudinal SST gradient in a positive feedback mechanism (Fedorov et al., 2010). In other words, source region SSTs are potentially relevant to tropical thermocline dynamics in terms of both stand-alone trends, and in the context of distant SST trends. Tropical cyclones force vertical mixing in the upper ocean, which deepens the tropical thermocline (Bueti et al., 2014; Jansen et al., 2010). Tropical cyclone activity was higher in the early Pliocene than at present, especially in the Pacific (Fedorov et al., 2010). Gradually increasing zonal and meridional SST gradients (Fedorov et al., 2015) would have promoted stronger Walker and Hadley circulation, thereby reducing tropical cyclone activity and raising the tropical thermocline (Brierley et al., 2009; Fedorov et al., 2010). The Pacific underwent a stronger reduction in tropical cyclone activity than the Atlantic (Fedorov et al., 2010), which may be linked to the contrasting evolution of the Pacific and Atlantic tropical thermocline depth. This feedback mechanism between midlatitude temperatures and tropical cyclone activity could further be explored in relationship to tropical thermocline depth. The hypothesized link should be tested in model simulations and could additionally be supported or refuted by better proxy records from midlatitudes.

5 Conclusions

Our new $\mathrm{Mg}/\mathrm{Ca}$, δ¹⁸O, and δ¹³C records from Site 959 indicate late Pliocene thermocline and nutricline deepening in the Eastern Equatorial Atlantic, starting at the end of the mPWP. A temporal discrepancy between Δδ¹⁸O and $\mathrm{\Delta }\mathrm{Mg}/\mathrm{Ca}$ can be explained by a transient reduction of the vertical salinity gradient in response to intensification of west African monsoonal precipitation during the mPWP. Tropical thermocline deepening occurred across the Atlantic basin, as is indicated by a nearly identical Δδ¹⁸O record at Site 1000 in the Caribbean Sea. This seems inconsistent with the hypothesis that Pliocene global cooling was linked to tropical thermocline shoaling, as observed in the Pacific, although Atlantic thermocline adjustments may have had less of an effect on the global heat budget because the Atlantic basin is smaller. The mechanisms driving Pliocene thermocline changes are currently speculative. In climate models, Central American Seaway closure leads to global thermocline shoaling, but during most of the early and late Pliocene, opposite changes occurred in the tropical Atlantic and Pacific thermoclines. We therefore suggest further exploring alternative mechanisms, such as source water temperature forcing in the Southern Ocean in combination with tropical cyclone activity. Due to basin morphology, tropical Atlantic thermocline waters are sourced from the Southern Hemisphere, whereas both hemispheres feed the tropical thermocline waters in the Pacific. Asymmetric Northern and Southern Hemisphere temperature trends may therefore have contributed to the differences in Atlantic and Pacific thermocline depth changes. A positive feedback loop between meridional SST gradients and tropical cyclones may have amplified vertical thermocline movements.