Plio-Pleistocene Perth Basin water temperatures and Leeuwin Current dynamics (Indian Ocean) derived from oxygen and clumped-isotope paleothermometry

The Pliocene sedimentary record provides a window into Earth's climate dynamics under warmer-than-present boundary conditions. However, the Pliocene cannot be considered a stable warm climate that constitutes a solid baseline for middle-of-the-road future climate projections. The increasing availability of time-continuous sedimentary archives (e.g., marine sediment cores) reveals complex temporal and spatial patterns of Pliocene ocean and climate variability on astronomical timescales. The Perth Basin is particularly interesting in that respect because it remains unclear if and how the Leeuwin Current sustained the comparably wet Pliocene climate in Western Australia, as well as how it influenced Southern Hemisphere paleoclimate variability. To constrain Leeuwin Current dynamics in time and space, this project obtained eight clumped-isotope Δ₄₇ paleotemperatures and constructed a new orbitally resolved planktonic foraminifera (Trilobatus sacculifer) stable isotope record (δ¹⁸O) for the Plio-Pleistocene (4–2 Ma) interval of International Ocean Discovery Program (IODP) Site U1459. These new data complement an existing TEX₈₆ record from the same site and similar planktonic isotope records from the Northern Carnarvon Basin (Ocean Drilling Program (ODP) Site 763 and IODP Site U1463). The comparison of TEX₈₆ and Δ₄₇ paleothermometers reveals that TEX₈₆ likely reflects sea surface temperatures (SSTs) with a seasonal warm bias (23.8–28.9 ^∘C), whereas T. sacculifer Δ₄₇ calcification temperatures probably echo mixed-layer temperatures at the studied Site U1459 (18.9–23.2 ^∘C). The isotopic δ¹⁸O gradient along a 19–29^∘ S latitudinal transect, between 3.9 and 2.2 Ma, displays large variability, ranging between 0.5 ‰ and 2.0 ‰. We use the latitudinal δ¹⁸O gradient as a proxy for Leeuwin Current strength, with an inverse relationship between both. The new results challenge the interpretation that suggested a tectonic event in the Indonesian Throughflow as the cause for the rapid steepening of the isotopic gradient (0.9 ‰ to 1.5 ‰) around 3.7 Ma. The tectonic interpretation appears obsolete as it is now clear that the 3.7 Ma steepening of the isotopic gradient is intermittent, with flat latitudinal gradients (∼0.5 ‰) restored in the latest Pliocene (2.9–2.6 Ma). Still, the new analysis affirms that a combination of astronomical forcing of wind patterns and eustatic sea level controlled Leeuwin Current intensity. On secular timescales, a period of relatively weak Leeuwin Current is observed between 3.7 and 3.1 Ma. Notably, this interval is marked by cooler conditions throughout the Southern Hemisphere. In conclusion, the intensity of the Leeuwin Current and the latitudinal position of the subtropical front are both long-range effects of the same forcing: heat transport through the Indonesian Throughflow (ITF) valve and its propagation through Indian Ocean poleward heat transport. The common ITF forcing explains the observed coherence of Southern Hemisphere ocean and climate records.

1 Introduction

Arid climates reign in large portions of the Australian continent, yet along the coastlines, the habitability of the Australian continent is often facilitated by moisture transport, associated with the many boundary currents surrounding the continent. The Leeuwin Current is one of those boundary currents, flowing southward along the Northwest Shelf of Australia (NWS), across the Carnarvon and Perth basins, ultimately rounding Cape Leeuwin to flow further east into the Great Australian Bight. Leeuwin Current intensity varies throughout the year, with strongest flow in austral winter when the latitudinal steric height gradient is steepest (Fig. 1) (Cresswell and Peterson, 2009; Waite et al., 2007; Godfrey and Ridgway, 1985; Ridgway and Godfrey, 2015; Church et al., 1989; Pearce, 1991). Thereby, the Leeuwin Current causes thermocline deepening, while also acting as a moisture and heat source for the Mediterranean-like climate around Perth, characterized by wet and mild winters. Throughout the Neogene, rainfall patterns in Western Australia experienced severe regime shifts on million-year timescales. Groeneveld et al. (2017) describe an extremely arid middle Miocene based on sabkha-like sediments found on the Northwest Shelf of Australia. Tagliaro et al. (2018) corroborate this climate interpretation by proposing sea level fall and regional aridity as causes for Miocene karst in the Northern Carnarvon Basin. Throughout the late Miocene, a northward shift of the westerlies allowed for a progressive increase in precipitation in southwest Australia (Groeneveld et al., 2017). At the same time, the siliciclastic Bare Formation was deposited on the NWS in a phase of deltaic margin progradation, thus correlating with increasing humidity in the northwest Australian hinterland (Tagliaro et al., 2018). By the latest Miocene and throughout most of the Pliocene, Western Australia was entirely governed by a wet or seasonally wet climate, i.e., the so-called “Humid Interval” (Christensen et al., 2017; Karatsolis et al., 2020). The switch from an arid Miocene to a wet Pliocene is also reflected in fossil pollen records from southern Australia (Sniderman et al., 2016). Finally, in the late Pliocene and early Pleistocene, climate transitioned back towards a more arid and more seasonal precipitation regime, albeit with important glacial–interglacial variability and more humid conditions during interglacials (Stuut et al., 2019; Fujioka et al., 2009; Stuut et al., 2014; Kuhnt et al., 2015; Gallagher et al., 2014; Dodson and Macphail, 2004; Dodson and Ramrath, 2001; He and Wang, 2021). The Late Pliocene aridification of Australia promoted C₄ over C₃ photosynthesis and probably constitutes the underlying reason for the late C₄ expansion in Australia, compared to other continents (Andrae et al., 2018).

Figure 1 Present-day oceanography of the eastern Indian Ocean. The southward flowing Leeuwin Current is strongest in late autumn and winter. During that season, the Leeuwin Current generates large eddies, causing significant deepening of the mixed layer (down to ∼200 m) and facilitating enhanced primary productivity (also see Fig. B1). During austral summer, the Leeuwin Current is weaker, with a corresponding thinner mixed layer and reduced productivity. (a) Sea surface temperature on 1 July 2021 (austral winter). (b) Sea surface temperature on 1 January 2022 (austral summer). (c) Time–depth potential temperature plot for the Site U1459 locality. Sea surface temperatures vary seasonally between 20 and 24 ^∘C. Figures were generated using EU Copernicus Marine Service Information with data from GLOBAL_ANALYSIS_FORECAST_PHY_001_024 (Zammit-Mangion and Wikle, 2020).

While the million-year-scale Neogene hydroclimate evolution of Western Australia is relatively well constrained and substantiated by multiple lines of evidence, important open questions remain regarding the temporal and causal relationships between Australian hydroclimate evolution and Leeuwin Current dynamics. Cane and Molnar (2001) hypothesize that east African aridification at 3–4 Ma was caused by Indonesian Throughflow (ITF) restriction, assuming a switch in the source of surface flow through the Indonesian seaway from South to North Pacific waters. This publication triggered an increased research interest in the causal relationship between Indo-Pacific paleoceanography and paleoclimate change. Karas et al. (2009) supports this hypothesis but points out that the switch likely took place at the subsurface level rather than at the surface. Using a general circulation model, Krebs et al. (2011) links the end of the “Humid Interval” and the observed Late Pliocene desertification of Australia to a reduction in ITF transport. In the same year, Karas et al. (2011a) were the first to explicitly link the Late Pliocene aridification of Western Australia to Leeuwin Current weakening in response to a reduced ITF. In 2015, International Ocean Discovery Program (IODP) Expedition 356 “Indonesian Throughflow” cored seven sites along Australia's western margin and obtained long, time-continuous sediment sequences that have been used to chart Pliocene-to-recent Leeuwin Current dynamics on astronomical timescales and to link those dynamics to the climate evolution of Western Australia. Auer et al. (2019) and De Vleeschouwer et al. (2018b, 2019b) focus on Late Pliocene glacial–interglacial variability around Marine Isotope Stage M2 (MIS M2, 3.3 Ma) and suggest that Pliocene aridity is punctuated during glacials when the Leeuwin Current was relatively weak. Contrary to this, during strong interglacials, winter aridity in northwest Australia is alleviated by a more active Leeuwin Current, which acts as a moisture source. Organic geochemistry TEX₈₆ studies reveal three important sea surface temperature (SST) cooling steps at sites under the influence of the Leeuwin Current at 3.3–3.1, 1.7–1.5, and 0.65 Ma (Smith et al., 2020; Petrick et al., 2019). These authors interpreted the cooling steps to result from a progressive weakening of the Leeuwin Current in response to ITF constriction, at least during glacial periods. This interpretation neatly accords with the gradual aridification of Western Australia during this period. Nevertheless, He et al. (2021) challenge this paradigm: they use zonal and meridional intersite temperature gradients and productivity reconstructions to postulate that the Leeuwin Current became stronger – not weaker – after the Mid-Pleistocene Transition. These authors interpret their TEX₈₆ records as a subsurface signal though, and further work is needed to test the He et al. (2021) interpretation. First, the increasing discrepancies between ${\mathrm{U}}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$ and TEX₈₆ at their studied Site U1461 over the last 1 Myr needs to be clarified, and the lack of correlation between their records and the benthic isotope stack needs to be explained.

On glacial–interglacial timescales, the Pleistocene behavior of the Leeuwin Current is comparable to that of the Pliocene, as it oscillates between two contrasting states. Interglacials were characterized by a stronger Leeuwin Current as higher eustatic sea levels allow for better connected shallow ITF pathways, for example, sourcing the Leeuwin Current over a flooded Sahul Shelf in northern Australia. Contrary, glacials were marked by shelf exposure along the coast of Western Australia, reduced ITF volume, and a stronger West Australian Current flowing opposite to the Leeuwin Current, all contributing to a weakening of the Leeuwin Current (Spooner et al., 2011; Wyrwoll et al., 2009; Petrick et al., 2019). Auer et al. (2021) refine this picture: a glacial sea level drop exceeding −45 m would expose the Sahul Shelf and cut off the Leeuwin Current from an important shallow-water source area. A weaker Leeuwin Current, in turn, would allow for the upwelling of sub-Antarctic Mode Water onto the Australian shelf and ultimately enhanced organic carbon burial. The coupling between Leeuwin Current intensity and glacial eustasy has been established by De Vleeschouwer et al. (2019b), using planktonic δ¹⁸O gradients along the Leeuwin Current pathway (19–29^∘ S) as an inverse proxy for the Leeuwin Current strength: a steep latitudinal gradient is observed when the Leeuwin Current is weak during glacials, and vice versa during interglacials (De Vleeschouwer et al., 2019b). However, the modeling results displayed in Fig. 2 suggest quite the opposite. When the Pliocene glacial simulation “Large M2” (Fig. 2b; Dolan et al., 2015) is compared to the warm-orbit Pliocene simulation “PlioMax” (Fig. 2a; Prescott et al., 2018), glacial cooling in the Carnarvon Basin seems to be more severe than in the Perth Basin (Fig. 2c). This pattern implies a flatter SST gradient during glacials compared to interglacials and is thus the exact opposite of what is observed in proxy data. In this case, it is likely that these model simulations misrepresent past ocean dynamics because they do not consider the effects of shelf exposure and a sea-level-induced reduction in shallow-water ITF transport. Indeed, relative sea level plays an important role in feeding and directing the Leeuwin Current, with proxy data evidencing that Leeuwin Current intensity is sensitive to small changes in boundary conditions.

Figure 2 Pliocene modeled SSTs. (a) Warm Pliocene simulation “PlioMax” for Marine Isotope Stage (MIS) K1 (3.06 Ma orbit and 405 ppm CO₂; Prescott et al., 2018). (b) Cold Pliocene simulation “Large M2” for MIS M2 (3.30 Ma orbit and 280 ppm CO₂; Dolan et al., 2015). (c) The SST difference between the cold and warm Pliocene simulations suggests a reduction of the latitudinal SST gradient along the Leeuwin Current (LC) pathway during glacials. This modeling result does not match proxy results. Paleoceanographic proxy data from De Vleeschouwer et al. (2019b) suggest a steeper SST gradient along the LC pathway during glacials compared to interglacials. This misrepresentation in the model likely occurs because these particular models do not consider the effects of shelf exposure and a sea-level-induced reduction in shallow-water ITF transport. Both processes are crucial in understanding climate and ocean dynamics on secular and astronomical timescales. NEC is the North Equatorial Current. SEC is the South Equatorial Current. WAC is the West Australian Current.

This paper evaluates the driving factors of Leeuwin Current intensity on astronomical timescales by extending the approach to calculate planktonic δ¹⁸O gradients along the Leeuwin Current pathway across the Plio-Pleistocene boundary. Moreover, it compares TEX₈₆ and clumped-isotope paleothermometers at IODP Site U1459 to constrain absolute temperatures (and eventually heat transport) throughout the Plio-Pleistocene. Both paleothermometers have their pitfalls: the origin of the TEX₈₆ signal is still debated, while the ¹³C¹⁸O¹⁶O isotopologue is so uncommon that clumped-isotope paleothermometry still deals with large uncertainties. Per contra, the uncertainty on clumped-isotope measurements is largely random, and therefore averaged measurements provide meaningful results, even when error bars are wide. This makes clumped-isotope paleotemperatures a useful benchmark for comparisons with proxies that potentially face systematic biases (e.g., TEX₈₆). By comparing different ocean temperature reconstructions, our ambition is to obtain a better understanding of the links between Leeuwin Current dynamics and the climate evolution of the Southern Hemisphere.

2 Materials and methods

2.1 IODP Site U1459 (Perth Basin)

IODP Site U1459 (192 m water depth) lies in the northern part of the Perth Basin and is the southernmost site drilled in the IODP Expedition 356 latitudinal transect (Figs. 1, 2). The main objective of drilling Site U1459 was to obtain insight into the pre-Quaternary history of the Houtman Abrolhos reef complex and Leeuwin Current. In this work, the portion of Site U1459 between 55.65 and 106.41 m composite depth (mcd) is studied, which corresponds to lithostratigraphic Units III and IV in Gallagher et al. (2017). Unit IV is an unlithified, cream to light brown packstone, ascribed to an outer shelf to upper bathyal area. Unit III is an alternatingly light gray – greenish-gray unlithified to partially lithified packstone to grainstone. Compared to Unit IV, Unit III is defined by a marked increase in glauconite content and bioclasts (benthic foraminifera, bivalves, echinoderms, bryozoans, and gastropods). Benthic foraminifera assemblages ascribe Unit III to a middle to outer shelf environment.

De Vleeschouwer et al. (2019b) published several proxy series for Site U1459, including X-ray fluorescence (XRF)-derived element ratios, TEX₈₆ paleotemperatures, and planktonic δ¹⁸O and δ¹³C. The cyclostratigraphic age model in that paper consists of 15 age–depth tie points between 64 and 161 mcd, corresponding to 2.51 and 5.46 Ma. Here, we extend the studied interval into the early Pleistocene, which necessitates the addition of three age–depth tie points using a similar tuning approach as in De Vleeschouwer et al. (2019b). As a result, the U1459 age–depth model now extends to 51.42 mcd, equivalent to 1.93 Ma (Table A1; see also Sect. 3.1).

During IODP Expedition 356, natural gamma radiation (NGR) spectra of the recovered sediments were measured using the whole-round section logger (Gallagher et al., 2017). These spectra can be used to quantify sedimentary uranium and thorium contents through a freely available algorithm published by De Vleeschouwer et al. (2017).

2.2 IODP Site U1463 and ODP Site 763 (Northern Carnarvon Basin)

All new proxy records presented in this manuscript come from IODP Site U1459. In addition, we use previously published proxy records from IODP Site U1463 and Ocean Drilling Program (ODP) Site 763 to calculate planktonic δ¹⁸O and SST gradients along the Australian West Coast. IODP Site U1463 (145 m water depth) was cored in the Northern Carnarvon Basin (Beagle sub-basin) in 2015 and quickly became a regional reference site for reconstructing Miocene-to-recent climate, ocean, and basin dynamics (e.g., Christensen et al., 2017; De Vleeschouwer et al., 2018b; Auer et al., 2019; Smith et al., 2020; Karatsolis et al., 2020; Mccaffrey et al., 2020; Tagliaro et al., 2018; Gurnis et al., 2020; Groeneveld et al., 2021). The Pliocene interval of U1463 consists of homogeneous fine-grained mudstones with subordinate wackestone and packstone intervals, with benthic foraminifera assemblages indicative of a middle to outer shelf environment. ODP Site 763 (1367 m water depth) was also cored in the Northern Carnarvon Basin (Exmouth Plateau) in 1988 and consists of light-gray to white foraminifer–nannofossil ooze (Haq et al., 1990). The age–depth model is from Karas et al. (2011a), with a minor adjustment in the youngest part of the record (Table A2).

Both sites are located within the trajectory of the Leeuwin Current, at 18.96 and 20.58^∘ S, respectively (∼595 km apart). Because of their proximity and congruous oceanographic setting, we combine proxy series from both sites to construct a “Northwest Shelf of Australia” end-member to compare the newly generated U1459 data with and to calculate paleoceanographic gradients along the Leeuwin Current pathway.

2.3 Oxygen and carbon stable isotope analyses

We present a new dataset of 238 stable (δ¹³C and δ¹⁸O) isotope analyses, measured on calcite tests of the shallow dwelling planktonic foraminifer Trilobatus sacculifer. The new dataset connects to the younger end of the previously published Pliocene isotope record in De Vleeschouwer et al. (2019b; N=143) and extends it into the Pleistocene. The new measurements come from the stratigraphic interval between 82.9 and 55.65 mcd, and sampling occurred at a median spatial resolution of 15 cm according to the shipboard composite depth (mcd) scale, which resulted from a color-reflectance-based correlation between holes U1459A and U1459B (Gallagher et al., 2017). Specimens were picked from the 315–355 µm-size fraction to avoid size effects in δ¹⁸O values (Elderfield et al., 2002) and to be methodologically consistent with other planktonic isotopic records in the study area (De Vleeschouwer et al., 2018b; 2019b; Karas et al., 2011a). During picking, specimens with gametogenic calcite final chambers and partly dissolved specimens were avoided. All samples were measured using a Finnigan MAT 252 gas isotope ratio mass spectrometer connected to a Kiel III automated carbonate preparation device at the Center for Marine Environmental Sciences (MARUM). Isotopic data are reported in standard delta notation vs. Vienna Pee Dee Belemnite (V-PDB). We calibrated all measurements against the in-house standard (homogenized Solnhofen limestone powder), which in turn is calibrated against the NBS-19 reference material. Over the measurement period, the standard deviations (1σ) of the in-house standard (N=90) were 0.04 ‰ for δ¹³C and 0.06 ‰ for δ¹⁸O.

2.4 Clumped-isotope thermometry in planktonic foraminifera

Clumped-isotope measurements were carried out using calcite tests of T. sacculifer in the 315–355 µm size fraction. Specimens from up to four adjacent samples were pooled to obtain sufficient material for clumped-isotope paleothermometry (150–350 specimens, 2–5 mg). Potential contaminants were then removed using a modified version of the cleaning protocol for foraminiferal $\mathrm{Mg}/\mathrm{Ca}$ analysis of Barker et al. (2003), without the H₂O₂ treatment steps (Grauel et al., 2013; Peral et al., 2018). The clumped-isotope analysis of the cleaned carbonate powders took place at the Vrije Universiteit Brussel, Belgium (AMGC-VUB lab), using a Nu Instruments Perspective-IS stable isotope ratio mass spectrometer (SIRMS) in conjunction with a Nu-Carb carbonate sample preparation system. This setup also includes a fully automated adsorptive trap purification system to remove contaminants. Between 500 and 600 µg of carbonate powder reacts for 10 min at 70 ^∘C after the automated injection of 120 µL of H₃PO₄. The CO₂ gas produced is expanded into a water trap for 5 min, held at −95 ^∘C. Subsequently, the CO₂ gas moves through a PoraPak™ Q-packed adsorption trap, held at −34 ^∘C for 25 min, and trapped into a liquid nitrogen (LN2) cold finger. After a yield measurement utilizing a pressure transducer, the sample is transferred into the sample cold finger, for 3 min, within the inlet system of the mass spectrometer. The cold finger on the sample side of the dual inlet acts as a constant depletion volume once the gas has been expanded. This is matched by an identical depletion volume on the reference side of the dual-inlet system. This allows for a constant depletion rate of both sample and reference gas during the data acquisition. Additionally, the absorptive trap is cleaned/degassed by active vacuum pumping at 150 ^∘C between sample extractions for 25 min to remove any contaminants that may have been trapped during the gas transfer. Reference and sample gases are alternatingly measured on six Faraday collectors ($m/z$ 44–49) and analyzed in three “blocks” of 20 cycles each with a counting time of 20 s. Measurements sum up to 20 min of integration time per replicate. A zero (no-gas) background measurement and automatic peak centering (on $m/z$ 45) are performed at the start of each sample measurement. The reference gas beam is pressure balanced to the sample gas beam and depletes evenly through matched length capillaries; initial beam balance ($m/z=\mathrm{44}$) is set to 80 nA, weakening to approximately 45 nA. The reference side of the dual inlet was refilled with reference gas every seven analyses. Total analysis time (including the reaction, PoraPak purification, and integration time) is approximately 1 h 20 min per sample. The subsequent sample starts its preparation while the mass spectrometer is still analyzing the previous sample. Possible contamination is monitored on ETH1–4 standards by scrutinizing Δ48 raw and Δ49 raw values for extraordinary deviations from the mean (Meckler et al., 2014). The ETH standards were measured following the recommendations of Kocken et al. (2019) with a sample to standard ratio of 1:1. Analyses and results are monitored in the lab using the Easotope software (John and Bowen, 2016). The raw Δ₄₇ values of sample-derived CO₂ were converted to the I-CDES 90 ^∘C scale, using the most recent values for the ETH-1, ETH-2, and ETH-3 carbonate reference materials (Bernasconi et al., 2021). The raw measured Δ₄₇ values were processed using the IUPAC isotopic parameters (Brand et al., 2010; Petersen et al., 2019; Daëron et al., 2016) within the ClumpyCrunch software (Daëron, 2021). The reported analytical uncertainties have been calculated using the “pooled” standardization method described by Daëron (2021), incorporating constraints on reproducibility available from both standard and sample analyses. A total of 307 Δ₄₇ replicate analyses (standards and samples) have been carried out in the framework of this study, of which one replicate analysis has been rejected as an outlier because it was >4σ away within entire multi-sample population, as well as within the population of PB06 replicates.

The T. sacculifer clumped-isotope measurements from Site U1459 were converted into calcification temperatures, using three different calibration schemes, which all follow the form of Eq. (1).

$$\begin{array}{}\text{(1)}& {\mathrm{\Delta }}_{\mathrm{47}}={\displaystyle \frac{\text{A}}{T^{\mathrm{2}}}}+\text{B}\end{array}$$

The most suitable calibration for the dataset presented here is that of Meinicke et al. (2020), as recalculated in Meinicke et al. (2021), because it is the most comprehensive foraminifera-specific calibration dataset currently available. Another relevant foraminifera-specific calibration for this study is the clumped-isotope paleothermometry calibration by Peral et al. (2018), as recalculated in Peral et al. (2022). In fact, the Meinicke et al. (2020, 2021) calibration fully includes the datasets of earlier calibrations by Peral et al. (2018, 2022) and Piasecki et al. (2019). The most recent Anderson et al. (2021) and Fiebig et al. (2021) calibrations have the advantage of covering a wide range of temperatures (0.5–1100 ^∘C). However, due to the ongoing debate regarding linear vs. polynomial behavior of the clumped-isotope calibration across a wide temperature range, we prefer working with the foraminifera-specific calibration that unites all data of Meinicke et al. (2020), Piasecki et al. (2019) and Peral et al. (2018). Therewith, we focus on the natural range of ocean temperatures, appropriate for the paleoceanographic purposes of this study. It is important to emphasize that the Peral et al. (2018) and Meinicke et al. (2020) regression parameters have been updated to match the new I-CDES 90 ^∘C scale (Bernasconi et al., 2021). The updated version of those calibrations can be found in Peral et al. (2022) and Meinicke et al. (2021), respectively, as well as in Table 1.

Table 1 Clumped-isotope paleothermometry calibration parameters for Eq. (1), using the I-CDES 90 ^∘C scale. The A and B parameters of the Meinicke et al. (2020) calibration were recalculated from the CDES to the I-CDES 90 ^∘C framework by Meinicke et al. (2021). The Peral et al. (2018) calibration parameters were recalculated by Peral et al. (2022). The reader is referred to Peral et al. (2022) for more information on the uncertainty on this calibration. A recalculation of the Anderson et al. (2021) calibration was not necessary, as they have been originally defined in the I-CDES 90 ^∘C framework. Reported uncertainties refer to ±1 SE.

Temperature uncertainties for individual samples in this study (68 % and 95 % confidence intervals) include both analytical and calibration uncertainty, and were calculated using a Monte Carlo approach.

3 Results

3.1 High-resolution oxygen isotope record and construction of the age–depth model

The new planktonic δ¹⁸O depth series (obtained from 55.65 to 82.90 mcd) smoothly connects to the previously published dataset (83.4–106.41 mcd; De Vleeschouwer et al., 2019b): the two datasets align in terms of absolute values at the transition point, and the variance and isotopic range also remain similar between the younger and older interval (Fig. 3a). The previously published dataset exhibits positive δ¹⁸O values up to 0.68 ‰ around 88 mcd, which were linked to the global MIS M2 glaciation by De Vleeschouwer et al. (2019b). A gradual trend towards more negative values characterizes the record in the few meters above this level. The new dataset starts at about the point where the gradual trend is concluded. The first ∼17 m of the new dataset (82.9–65.5 mcd) are marked by relatively negative (between 0 and −0.5 ‰) and low-variance isotope values. No distinct positive excursions occur in this interval. At 65.5 mcd, a rapid shift to more positive δ¹⁸O values occurs and this shift is accompanied by increased variability. This stratigraphic interval lies outside the range of the age–depth model that was constructed in De Vleeschouwer et al. (2019b). Hence, to convert depth to age, it was necessary to extend the available age–depth model by three additional age–depth tie points. It was not possible to use the $\mathrm{Ca}/\mathrm{Fe}$ depth series for this purpose because it has large data gaps in the upper interval (Fig. 3b). Instead, we aligned the planktonic δ¹⁸O record directly with a Southern Hemisphere eccentricity-tilt-precession (ETP) composite (Laskar et al., 2004) and the LR04 benthic stack (Lisiecki and Raymo, 2005). As we look at astronomical forcing signatures within the Site U1459 δ¹⁸O signal and in the latitudinal isotopic gradient along the western coast of Australia later in this work, we recognize that this age–depth modeling approach carries some risk for circular reasoning. However, no viable alternative technique allows the construction of a reliable age–depth model with 10⁵-year temporal resolution (i.e., an order of magnitude greater than the biostratigraphic time resolution in this interval). Nevertheless, we minimize the potential impact of circular reasoning by limiting the alignment of the isotope record to the target curves to only two tie points (one at MIS 92 and one at MIS 87). The third tie point is a biostratigraphic tie point at 1.93 Ma and thus does not impose any circular reasoning (Fig. 3d).

Figure 3 IODP Site U1459 depth series and age–depth model. (a) Planktonic δ¹⁸O_calcite depth series, produced by merging the measurement series in De Vleeschouwer et al. (2019b, N=143) with the new data presented here (N=238). The stratigraphic positions of marked glacial intervals (MIS 96, 100, M2, Gi-2, and Gi-4) are indicated. (b, c) XRF-derived $\mathrm{Ca}/\mathrm{Fe}$ and $\mathrm{K}/\mathrm{Al}$ depth series, interpreted in De Vleeschouwer et al. (2019b) as proxies for detrital input and aeolian kaolinite flux, respectively. Both element ratios are plotted on a logarithmic scale. (d) The age–depth model by De Vleeschouwer et al. (2019b) is extended into the early Pleistocene by three additional age–depth tie points. Biostratigraphic markers are indicated by brown diamonds. Age–depth model tie points are indicated by open squares.

3.2 Clumped-isotope-based calcification temperature reconstructions

All clumped-isotope measurements (standards and samples) carried out in the framework of this study are summarized in Table 2. The ETH-1, ETH-2, and ETH-3 standards define the I-CDES 90 ^∘C framework and, therefore, correspond to the values fixed by the InterCarb initiative (Bernasconi et al., 2021). The ETH-4 standard is used for quality control, as its Δ₄₇ composition is known, but it is not used to define the I-CDES 90 ^∘C isotopic framework. The measured value for the ETH-4 standard matches the expected value within error margins (±1 SE; Table 2). Throughout the studied interval of IODP Site U1459, eight samples have been measured for clumped-isotope paleothermometry (PB01 through PB08). Sample measurements resulted in Δ₄₇ compositions that range between 0.6018 ‰ and 0.6156 ‰. The availability of sample material directly determined the number of repeated aliquot measurements carried out per sample. Sample availability remains the main limiting factor for clumped-isotope paleothermometry, as ∼500 µg of sample material is required per aliquot in the AMGC-VUB lab.

Table 2 Clumped-isotope results. ETH-1, 2, and 3 define the I-CDES 90 ^∘C framework. Hence their Δ₄₇ are identical to the values reported by Bernasconi et al. (2021). The ETH-4 standard is used for monitoring measurement accuracy. The availability of sample material allowed for 6 to 12 repeated measurements (∼500 µg each). The standard error on sample measurements is calculated using the pooled standardization method described by Daëron (2021). The δ¹⁸O and δ¹³C sample values reported in this table were measured simultaneously with Δ₄₇ in the AMGC-VUB lab.

Table 3 Clumped-isotope paleothermometry. Clumped isotopes reconstructed calcification temperatures according to three different calibrations. The Meinicke et al. (2020) calibration is adopted throughout this work. The δ¹⁸O_sw of seawater was reconstructed by solving the Erez and Luz (1983) equation, using the clumped-isotope calcification temperature and the foraminiferal δ¹⁸O. Note that all uncertainties in this table refer to 95 % confidence intervals (CIs): in the case of calculated clumped-isotope paleotemperatures, these CIs reflect Monte Carlo propagated analytical and calibration uncertainties.

Figure 4 (a) Ocean temperature reconstructions based on TEX₈₆ (De Vleeschouwer et al., 2019b) and clumped isotopes (this study) are offset by about 3 ^∘C. Clumped-isotope error bars represent ±68 % (bold lines) and 95 % confidence intervals (thin lines). (b) Comparison of U1459 clumped-isotope calcification temperatures and the planktic δ¹⁸O_calcite time series. While clumped-isotope measurements quantify past ocean temperatures, the δ¹⁸O_calcite time series allows for an investigation of past oceanographic variability on secular and astronomical timescales.

The Δ₄₇ clumped-isotope compositions translate to Perth Basin calcification paleotemperatures between 18.6 and 23.9 ^∘C when using the Meinicke et al. (2020) or Peral et al. (2018) calibrations. The Anderson et al. (2021) calibration produces paleotemperatures that are ∼0.7 to 1.6 ^∘C cooler (Table 3). Compared to the available TEX₈₆-based SST reconstruction for IODP Site U1459, the clumped-isotope temperatures are about 3 ^∘C cooler. This result is independent of which TEX₈₆ surface calibration is applied: BAYSPAR (Tierney and Tingley, 2014) or TEX${}_{\mathrm{86}}^{\mathrm{H}}$ (Kim et al., 2010) (Fig. 4a). The exact origin of the TEX₈₆ signal remains controversial though: it could be a surface, a subsurface, or a combined signal. In the ITF and Leeuwin Current study region, most authors used TEX₈₆ as a proxy for surface temperatures (Petrick et al., 2019; Smith et al., 2020; De Vleeschouwer et al., 2019b). However, He et al. (2021) used it as a proxy for subsurface temperatures and Meinicke et al. (2021) inferred a mixed surface–subsurface signal. In this paper, we start from the assumption that TEX₈₆ is a proxy for surface temperatures, but we add that the BAYSPAR subsurface calibration yields temperatures that are ∼3.3 ^∘C cooler compared to the surface calibration, while the Ho and Laepple (2016) subsurface calibration yields cool temperatures (13.7–15.9 ^∘C) that correspond to modern-day lower thermocline temperatures. All clumped-isotope derived calcification temperatures intersect with the 5 % confidence level of the TEX₈₆ BAYSPAR surface calibration (Fig. 4a). The clumped-isotope paleotemperatures are also compared to a foraminiferal δ¹⁸O-based SST reconstruction (Fig. 5b), using the Perth Basin seawater δ¹⁸O reconstruction presented in Fig. 5a as an assumption. The δ¹⁸O_sw has been obtained by superimposing the Plio-Pleistocene glacial contribution to δ¹⁸O_sw (Rohling et al., 2021) to the present-day isotopic composition of Perth Basin seawater (+0.8 ‰ Vienna Standard Mean Ocean Water (V-SMOW), dashed line in Fig. 5a). We assume local δ¹⁸O_sw changes, e.g., in response to salinity changes, to be negligible in this context. This approach results in a high-resolution SST reconstruction with temperatures ranging between 16 and 24 ^∘C. Therewith, the δ¹⁸O-based reconstruction aligns with the clumped-isotope temperatures and the δ¹⁸O-based temperature reconstruction intersects with the ±1 SE intervals of all eight clumped temperatures (Fig. 5b).

Figure 5 (a) Seawater δ¹⁸O_sw reconstruction for the Perth Basin, constructed by superimposing the Rohling et al. (2021) glacial contribution to δ¹⁸O_sw on the present-day isotopic composition (+0.8 ‰ V-SMOW). (b) The U1459 planktic δ¹⁸O time series can be converted to a high-resolution calcification temperature time series by using the Erez and Luz (1983) equation and adopting the δ¹⁸O_sw reconstruction in panel (a). The high-resolution δ¹⁸O-derived temperature series is compared to the lower-resolution Δ₄₇-derived paleotemperatures. Clumped-isotope error bars represent ±68 % (bold lines) and 95 % confidence intervals (thin lines). (c) Δ₄₇-derived paleotemperatures are compared to the U1459 core-based $\mathrm{U}/\mathrm{Th}$ series (note the reversed y axis). We use $\mathrm{U}/\mathrm{Th}$ as a proxy for C_org burial and upwelling of cold sub-Antarctic Mode Water (SAMW) under a weak Leeuwin Current regime (proxy justification in Auer et al., 2021). Enhanced upwelling at the Plio-Pleistocene boundary does not however seem to have a significant effect on clumped-isotope calcification temperatures.

4 Discussion

4.1 Discrepancy between TEX₈₆ sea surface and clumped-isotope calcification temperatures

The TEX₈₆ SST reconstruction at Site U1459 is a few degrees warmer than the simulated annual average sea surface temperatures of the Pliocene climate models shown in Fig. 2a–b (20–24 ^∘C). At the same time, the TEX₈₆ surface temperatures are clearly and significantly warmer than the clumped-isotope calcification temperatures obtained by analyzing planktonic foraminifer T. sacculifer. The ∼3 ^∘C temperature difference between these proxies demands a discussion of the processes that may underpin the discrepancy. The TEX₈₆ temperature reconstruction was assessed as reliable in De Vleeschouwer et al. (2019b) because glycerol dialkyl glycerol tetraethers (GDGTs) were confirmed to be of marine origin, and all control parameters had values within the recommended range for reliable paleo-SST reconstructions based on TEX₈₆. In this study, we adhere to this assessment. Indeed, the TEX₈₆ reconstruction for the Plio-Pleistocene with temperatures between 23.8 and 28.9 ^∘C seems reasonable compared to the present-day mean annual temperature of about 22.8 ^∘C (Reynolds et al., 2002). Nonetheless, we point out at least two potential warm biases in the U1459 TEX₈₆ SST reconstruction. First, TEX₈₆ SST reconstructions along the Leeuwin Current pathway might suffer from a seasonal bias towards summer temperatures, although the effect of seasonality on TEX₈₆ is poorly understood. A seasonal warm bias at Site U1459 is supported by the observation that core-top TEX₈₆ SST are 3 ^∘C warmer than modern mean annual SST at Site U1463 in the Carnarvon Basin, upstream the Leeuwin Current (Smith et al., 2020). A similar warm bias has been observed at the nearby Perth Basin Site U1460 by Benjamin Petrick (personal communication, 2019). The TEX₈₆ seasonal warm bias along the Leeuwin Current pathway could be linked to the oligotrophic conditions that characterize surface waters off the Australian west coast during summer (when Leeuwin Current is weak and the eddy-induced nutrient flux is minimum), allowing Thaumarchaeota to dominate and thrive (Guo et al., 2021). A second potential TEX₈₆ warm bias in the Perth Basin arises from the southward displacement of sinking particles, as Benthien and Müller (2000) proposed for the western Argentine Basin. Using a typical sinking rate of 100 m d⁻¹, a Leeuwin Current velocity of 50 cm s⁻¹, and a paleo-water depth of ∼500 m, we estimate the order of magnitude of possible southward displacement of particles at several hundred kilometers. This potential warm bias might not be negligible because of the rather steep latitudinal temperature gradient offshore southwestern Australia but is insufficient to explain temperature differences >1 ^∘C. Moreover, the T. sacculifer clumped-isotope paleothermometer would be subject to a similar, yet smaller, effect (Takahashi and Be, 1984). Hence, we conclude that, at Site U1459, TEX₈₆ is likely to suffer from a warm bias, mainly driven by a seasonal effect.

The clumped-isotope paleothermometer is also subject to potential overprinting and biases. Indeed, we cannot fully exclude a diagenetic impact on Δ₄₇ and δ¹⁸O. Dolomite is recognized throughout Site U1459 and becomes a stronger influence downcore. However, we do not observe a systematic change towards more positive δ¹⁸O values with increasing core depth. On the contrary, the δ¹⁸O data behave opposite to the expected diagenetic imprint in high-carbonate settings (e.g., Stainbank et al., 2020). Another possible cold bias resides with seafloor recrystallization. While partly dissolved specimens were dismissed during picking and a cleaning protocol has been in place to remove secondary calcite prior to Δ₄₇ analysis, this factor cannot be fully excluded. Nevertheless, our confidence in the Δ₄₇ and δ¹⁸O data is strengthened by the preservation of the selected T. sacculifer specimens (Fig. C1), which was classified as moderate to good (Gallagher et al., 2017), the internal coherence of the δ¹⁸O and Δ₄₇ series (Figs. 4b, 5b), and the observation that the calculated δ¹⁸O_sw values (Table 3) are within the expected range for the Pliocene Perth Basin (Fig. 5a). The isotopic consonance between δ¹⁸O and Δ₄₇ at Site U1459 constitutes an important indicator, suggesting that it is unlikely that isotopic scrambling (i.e., a reset of ¹³C–¹⁸O clumping) within the calcite of the analyzed foraminiferal tests occurred. A factor that might have played a role in causing the TEX₈₆ – clumped-isotope temperature offset is a difference in determinative water depth: clumped-isotope calcification temperatures likely reflect mixed-layer water temperatures rather than true surface temperatures. The mixed layer offshore southwest Australia can be fairly deep, down to 200 m during austral winter (Fig. 1c). This assessment is underpinned by results presented by Zhang et al. (2019). In a study in the Indonesian Throughflow region, they found that T. sacculifer calcified at the bottom of the mixed layer, or even at the upper thermocline. These authors suggest that the temperature difference between T. sacculifer and Globigerinoides ruber sensu stricto, a true surface-water dweller, could be used as a proxy for mixed-layer depth in the region. Their results were largely in line with earlier work by Rippert et al. (2016), who derived a T. sacculifer calcification depth of ∼120 m in the West Pacific Warm Pool, which is at the bottom of the mixed layer or the top of the thermocline. Nevertheless, the extension of the vertical habitat of T. sacculifer into the upper thermocline remains controversial, for example, because it contradicts an earlier study off Indonesia: Mohtadi et al. (2011) calculated much shallower calcification depths (between 20 and 75 m) for T. sacculifer within the mixed layer. Therewith, the Mohtadi et al. (2011) results are in agreement with a recent clumped-isotope study by Meinicke et al. (2021), who determined the apparent calcification depth of T. sacculifer in the Western Pacific Warm Pool to be around 75 m water depth, within the lower mixed layer. Regardless of the controversy, it becomes increasingly clear that T. sacculifer has a wider vertical range than G. ruber. For example, a large-scale analysis of plankton net haul data from the subtropical eastern North Atlantic reports average living depths for T. sacculifer between 15 and 200 m (Rebotim et al., 2017), and a study in the Kuroshio Current system also found T. sacculifer down to 200 m water depth (Kuroyanagi and Kawahata, 2004). Overall, these recent planktonic foraminiferal habitat studies support the notion that T. sacculifer is a species able to adapt its temporal (lunar cycle, seasonal cycle) and spatial (water depth) habitat to local oceanographic conditions (Jonkers and Kuèera, 2015, 2017; Kretschmer et al., 2018). In addition, isotopic proxies measured on T. sacculifer calcite tests could also be affected by a cold bias through the addition of gametogenic calcite. This effect has been described by Bé (1980), who noticed that living T. sacculifer shells collected during towing were smaller than surface sediment specimens. The size difference was ascribed by Duplessy et al. (1981) to the addition of calcite (gametogenic calcification) during their descent below the euphotic zone. This crust of secondary calcite is typically formed below the thermocline and can contribute up to ∼20 % of the post-gametogenic shell. Hence, gametogenic calcification can impose a significant cold bias on the isotopically reconstructed temperatures and result in unexpectedly deep apparent calcification temperatures (Wycech et al., 2018). However, we consider this effect to be limited thanks to the stringent selection of specimens during picking and the cleaning protocol.

We thus explain the temperature discrepancy between TEX₈₆ and clumped isotopes by a combination of two main processes: (1) T. sacculifer descending in the water column during austral winter, when the mixed layer is up to 200 m thick (Fig. 1), and (2) TEX₈₆, reflecting sea surface temperatures with a bias towards warmer summer temperatures.

In Fig. 5c, we compare the clumped-isotope temperatures to $\mathrm{U}/\mathrm{Th}$ ratios. Here, we employ the $\mathrm{U}/\mathrm{Th}$ ratio as a proxy for organic carbon burial, following the proxy interpretation by Auer et al. (2021) for the mid-Pleistocene interval of nearby Site U1460. Enhanced C_org burial during intense glacial stages is linked to the upwelling of nutrient-rich sub-Antarctic Mode Water onto the Western Australian Shelf. The upwelling of those waters triggered enhanced productivity and caused a contraction of the mixed layer and a marked steepening of the thermocline. With that perspective, it is interesting to note that glacial stages MIS Gi6, Gi4, G4, G2, 100, 98, 94, 92, and 88 are marked by high $\mathrm{U}/\mathrm{Th}$ (Fig. 5c). However, the upwelling of sub-Antarctic Mode Water at times of weak Leeuwin Current does not seem to have a strong effect on the clumped-isotope and δ¹⁸O-derived calcification temperatures for these glacials (Fig. 5b–c). This result indicates that T. sacculifer limited its habitat to the mixed layer, avoiding the thermocline at Site U1459 throughout the studied interval.

4.2 Isotopic gradients along the Leeuwin Current pathway as a measure for poleward heat transport

An effective way to investigate the temporal and spatial dynamics of the Leeuwin Current is by calculating proxy gradients at several localities along its pathway. In a previous paper, De Vleeschouwer et al. (2019b) utilized the δ¹⁸O gradient between the Northern Carnarvon Basin (sites U1463 and 763) and the Perth Basin (Site U1459) to reconstruct mid-to-late Pliocene (3.9–3.1 Ma) Leeuwin Current strength. He et al. (2021) adopted an analogous approach, calculating differences in TEX₈₆-based temperatures between Site U1461 (Northern Carnarvon Basin) and Site U1459 (Perth Basin). The low-resolution reconstruction of Leeuwin Current intensity by He et al. (2021) is in good agreement with our high-resolution δ¹⁸O gradient reconstruction (Fig. D1). Both proxies are completely independent of each other, implying that the covariation of their gradients is temperature driven and constitutes a proxy for Leeuwin Current strength.

Figure 6 (a) ODP Site 763 and IODP Site U1463 δ¹⁸O time series from the Northern Carnarvon Basin (purple), compared to IODP Site U1459 δ¹⁸O time series from the Perth Basin (orange). (b) Isotopic gradient between the Northern Carnarvon Basin (19^∘ S) and the Perth Basin (29^∘ S) along the Leeuwin Current pathway. (c–d) Forcing factors that are thought to influence Leeuwin Current intensity include eustatic sea level (dark green; Rohling et al., 2021; using the data in Westerhold et al., 2020) and the summer (21 June) intertropical insolation gradient (SITIG). They are combined with equal weights in a composite record after normalizing. (e–f) Multi-taper method power spectrum (Thomson, 1982) of the isotopic gradient and forcing composite. Their spectral signature is similar, dominated by obliquity and precession, and without a strong 100 kyr imprint. Spectral peaks in panel (e) are not strongly expressed because the underlying records have been subject to conservative astronomical tuning and age–depth modeling.

This paper extends the isotopic gradient approach across the Plio-Pleistocene boundary to 2.2 Ma (Fig. 6a, b). The extended analysis sheds new light on previous interpretations because, in the earlier work, a rapid increase in isotopic gradient (from 0.9 ‰ to 1.5 ‰, around 3.7 Ma) was interpreted as a permanent weakening of the Leeuwin Current. It was claimed that this rapid change was likely caused by a tectonic reorganization of the Indonesian Throughflow. However, the extension of the isotopic gradient analysis presented here reveals that the steepening of the isotopic gradient was not permanent. Indeed, one observes a gradual return to flat isotopic gradients throughout the Late Pliocene (3.1–2.9 Ma, Fig. 6b). By the latest Pliocene, between 2.9 and 2.6 Ma, latitudinal δ¹⁸O gradients were again as flat as 0.5 ‰ (Fig. 6b). Such flat gradient is indicative of a strong Leeuwin Current and was previously observed during the Early Pliocene, between 3.9 and 3.7 Ma. In other words, the rapid increase in isotopic gradient at 3.7 Ma and weakening of the Leeuwin Current were transient and are thus unlikely to have been caused by a tectonic reorganization in the throughflow region.

Figure 7 (a) Comparison of the isotopic gradient (dark gray) and forcing composite (green) reveals covariation on astronomical timescales. On secular timescales, a period of steep isotopic gradient and thus weak Leeuwin Current is observed between 3.7 and 3.1 Ma (gray rectangle) (b) Comparison of clumped-isotope calcification temperatures (blue crosses) and forcing composite (dark gray). The absence of temperature change between different samples indicates that the Leeuwin Current remained active all throughout the studied interval, despite variability in Leeuwin Current intensity.

The new analysis corroborates other interpretations proposed in De Vleeschouwer et al. (2019b): they argue that eustatic sea level and astronomical insolation forcing of wind patterns regulate Leeuwin Current dynamics on astronomical timescales. These two forcing factors are represented in Fig. 6c by token of the eustatic sea level reconstruction by Rohling et al. (2021) and the intertropical insolation gradient on 21 June (SITIG; Bosmans et al., 2015). Eustasy plays a role because sea level lowstands weaken ITF connectivity and reduce the availability of shallow Leeuwin Current source waters. Wind patterns are important because the Leeuwin Current has to overcome strong southerly winds during austral summer but benefits from southeasterly winds blowing off the Australian continent during winter. When the SITIG index is maximum, the Hadley cell over Australia is enhanced during austral winter but weakened during summer. Both effects spur Leeuwin Current flow under maximum SITIG. The sea level and insolation series in Fig. 6c have been normalized and added with equal weights to arrive at the composite signal in Fig. 6d. We thus propose a causal relationship between the combined forcing (Fig. 6d) and the isotopic gradient (Fig. 6b). This interpretation conforms with the broadly similar spectral character of both curves (Fig. 6e–f): both spectra are dominated by obliquity and precession and do not exhibit a 100 kyr eccentricity peak. We emphasize that the obliquity and precession spectral peaks in Fig. 6e occur at the expected frequencies, even though the underlying age–depth models (sites 763, U1463, and U1459) have been constructed using conservative tuning approaches. Indeed, none of those sequences have been tuned at the level of individual precession or obliquity cycles, with individual age–depth models not exceeding one tie point every ∼100 kyr. The most compelling argument for the proposed causal connection comes from the remarkable covariation between forcing and isotopic gradient on astronomical timescales, displayed in Fig. 7. Different glacial stages are recognizable in the forcing composite series and the isotopic gradient: intervals with a low forcing are generally characterized by a steep isotopic gradient and thus a weak Leeuwin Current. Also, the relative strength of individual glacial–interglacial cycle is consistent between forcing composite and isotopic gradient.

Figure 8 (a) Map of the South Atlantic, Indian, and Southern oceans indicating the location of DSDP, ODP, and IODP sites discussed in this work. LC is the Leeuwin Current; SEC is the South Equatorial Current; AC is the Agulhas Current; ARC is the Agulhas Return Current. (b–d) Different alkenone-based SST reconstructions (blue) in the Agulhas/Benguela region (Site 1087 from Petrick et al., 2018; Site 1090 from Martínez-Garcia et al., 2010; Site 1088 from Herbert et al., 2016) indicate late Pliocene cooling (gray rectangle), resulting from a more northerly position of the subtropical front (STF).

On longer secular timescales, the two curves seem to be decoupled in certain parts of the studied interval. Especially, between 3.7 and 3.1 Ma, the Leeuwin Current isotopic gradient is steeper than what one would derive from the similarity between forcing and gradient in the rest of the record (Fig. 7), with relatively low amplitude changes. In fact, the 3.7–3.1 Ma interval is characterized by isotopic gradients that are equally steep as the early Pleistocene interval between 2.6 and 2.2 Ma. Based on these results, we infer a prolonged period of Leeuwin Current weakness during the late Pliocene (3.7–3.1 Ma). Its timing chiefly corresponds to a time of Southern Hemisphere cooling, marked by the oldest tills in Argentinian Patagonia (3.6 Ma; Clague et al., 2020), a period of continental aridity in the south-central Andes (3.6–3.3 Ma; Amidon et al., 2017), the absence of the Dictyocha silicoflagellate and low abundance of calcareous nannoplankton between 3.5 and 3.2 Ma at ODP Site 751 on the Kerguelen Plateau (Bohaty and Harwood, 1998), a cooling of more than 3 ^∘C between 3.7 and 3.5 Ma based on warm vs. cold silicoflagellates from ODP Site 1165 on the Prydz Bay continental rise (Escutia et al., 2009), a marked decrease in sub-Antarctic zone diatoms because of surface ocean water cooling at Deep Sea Drilling Project (DSDP) Site 513 in the Atlantic sector of the Southern Ocean (3.6 Ma; Kato, 2020), and a gradual increasing trend in seasonal sea ice indicating diatoms at ODP Site 689 in the Weddell Sea (starting at 3.7 Ma; Kato, 2020) (see also Fig. 8).

The temporal evolution of the Leeuwin Current isotopic gradient also resembles alkenone-based SST reconstructions from the wider Agulhas region (South Atlantic, Fig. 8), which all show a general cooling trend between 3.7 and 3.1 Ma, followed by mostly stable and relatively warm temperatures between 3.1 Ma and the earliest Pleistocene (Site 1087 from Petrick et al., 2018; Site 1090 from Martínez-Garcia et al., 2010; Site 1088 from Herbert et al., 2016). The SST evolutions of those sites are coupled in varying degrees to latitudinal shifts in sub-Antarctic climate belts. The synchronicity between Leeuwin Current intensity and Agulhas SSTs thus suggests that there is a mechanistic link between the latitudinal position of the subpolar frontal system and Leeuwin Current intensity. This also implies that the rapid steepening of isotopic gradients at ∼3.7 Ma was probably associated with an intensified SAMW production in the Indian Ocean and possibly an equatorward shift of Southern Hemisphere climate belts rather than a tectonic event in the ITF region. Despite the covariation, open questions remain as to the exact cause and effect. It is rather unlikely that Leeuwin Current weakening would be the direct effect of an equatorward shift of the subtropical front. That is because such change would steepen the annual mean steric height gradient alongshore Western Australia, which would result in strengthening, not weakening, of the Leeuwin Current. The other way around, Karas et al. (2011a) suggested that a dwindling surface ITF has important implications for the Leeuwin Current and the poleward heat flux, possibly leading to a cooling of the Benguela upwelling system. Building on that suggestion, De Vleeschouwer et al. (2018b) described a positive feedback mechanism for Pliocene Southern Hemisphere climate change: an initial reduction in ITF heat transport leads to (sub-)Antarctic cooling, further glacio-eustatic sea level fall, and more Leeuwin Current weakening. Based on these hypotheses, it can be speculated that a reduction in poleward heat transport through the Leeuwin Current and the Agulhas (Return) Current caused the observed Southern Hemisphere cooling between 3.7 and 3.1 Ma. However, in this scenario, the question remains as to what triggered the strong cutback in poleward heat transport at ∼3.7 Ma in the first place. Possibly, the relatively large sea level drop related to the Early/Late Pliocene expansion of Northern Hemisphere ice sheets (De Schepper et al., 2014) was the trigger for the reduction in Indian Ocean (Leeuwin and Agulhas current) poleward heat transport. This initial change could have been amplified through the positive feedback mechanism in which the reduction of Indian Ocean poleward heat transport advances the thermal isolation of Antarctica (Karas et al., 2011a; De Vleeschouwer et al., 2018b), bringing the Southern Hemisphere into a long-lasting “cool state”.

On the other hand, the absolute TEX₈₆ and clumped-isotope paleotemperatures demonstrate that the Leeuwin Current continued to operate throughout the Pliocene, even when the latitudinal isotopic gradient became steeper. These new datasets thus provide important constraints for further ocean and climate modeling that can test the remaining open questions: were Pliocene sea level eustatic lowstands capable of sufficiently reducing Indian Ocean poleward heat transport so that it could influence Southern Hemisphere climate (without completely halting the Leeuwin Current)? Or were reductions in Leeuwin Current intensity the (indirect) result of an equatorward expansion of the sub-Antarctic climate belts, despite the corresponding steepening of the alongshore steric height gradient?

5 Conclusions

We have reconstructed a Plio-Pleistocene history of Leeuwin Current dynamics and Perth Basin ocean water temperatures, using a multiproxy approach at IODP Site U1459. Clumped-isotope and TEX₈₆-derived paleotemperatures for Site U1459 are about 3 ^∘C apart, which we explain by a combination of two processes: a seasonal warm bias on TEX₈₆ and a difference in the exact water depth represented by both proxies. Clumped isotopes measured on T. sacculifer probably mirror mixed-layer calcification temperatures, whereas TEX₈₆ is assumed to reflect surface waters.

The Plio-Pleistocene evolution of Leeuwin Current strength is assessed by calculating δ¹⁸O differences along a 19–29^∘ S latitudinal transect. Steep gradients are adopted as a proxy for weak Leeuwin Current intensities, and vice versa. This analysis confirmed the dependency of Leeuwin Current intensity on eustatic sea level and insolation forcing of wind patterns. Moreover, it indicates a “weaker than expected” Leeuwin Current between 3.7 and 3.1 Ma. This time interval chiefly corresponds to cool climate conditions across the Southern Hemisphere and a more northerly position of the subtropical front. The Leeuwin Current isotopic gradient exhibits remarkable congruity with SST records from the Southern Atlantic Ocean, supporting earlier proposed links between Leeuwin Current intensity and oceanographic change in the Agulhas/Benguela region. While this result underlines the importance of Indian Ocean poleward heat transport through the Leeuwin Current and the Agulhas (Return) Current, it remains an open question why the Leeuwin Current remained “weaker than expected” for ∼600 kyr after 3.7 Ma.

Appendix A: Age–depth models of IODP Site U1459 and ODP Site 763

Table A1 Age–depth model for the Plio-Pleistocene interval of IODP Site U1459.

Table A2 A minor chronology adjustment of Site 763 compared to Karas et al. (2011a) was needed in the youngest part of the record. This is because in the interval younger than 3.1 Ma, the Karas et al. (2011a) age–depth model is solely based on magnetic reversals. The minor chronology adjustment shown here was necessary to make sure that heavy planktonic isotope compositions correspond to glacials, and vice versa.

Appendix B: Present-day oceanography of the southeastern Indian Ocean

Figure B1 Present-day oceanography of the southeastern Indian Ocean. The southward flowing Leeuwin Current is strongest in late autumn and winter. During that season, the Leeuwin Current forms large eddies, causing significant deepening of the mixed layer (up to ∼200 m) and facilitating enhanced primary productivity (shown here are data from 1 July 2021). During austral summer, the Leeuwin Current is weaker (shown here are data from 1 January 2021), with a thinner mixed layer and reduced productivity is reduced. Figure was generated using EU Copernicus Marine Service Information with https://resources.marine.copernicus.eu/product-detail/GLOBAL_ANALYSIS_FORECAST_PHY_001_024/INFORMATION (last access: 22 May 2022) (Zammit-Mangion and Wikle, 2020) and https://resources.marine.copernicus.eu/product-detail/GLOBAL_ANALYSIS_FORECAST_BIO_001_028/INFORMATION (last access: 22 May 2022).

Appendix C: T. sacculifer scanning electron microscope images

Figure C1 Representative SEM images of T. sacculifer in sample U1459A 9FCC 9–24 cm, illustrating their preservation (Gallagher et al., 2017).

Appendix D: TEX₈₆ and isotopic gradients along the Leeuwin Current pathway as a proxy for Leeuwin Current intensity

Figure D1 Comparison of the TEX86 gradient between sites U1461 and U1459 (as in Fig. 3c in He et al., 2021), and the isotopic gradient from this study. Both proxies are completely independent of each other, yet show similar patterns throughout the Plio-Pleistocene, corroborating the assumption that TEX86 and isotopic gradients along the Leeuwin Current pathway are temperature driven and thus can serve as a proxy for Leeuwin Current intensity.