A Ternary Mixing Model Approach Using Benthic Foraminifer δ13C‐δ18O Data to Reconstruct Late Pliocene Deep Atlantic Water Mass Mixing

Late Pliocene deep Atlantic δ13C data have been interpreted as evidence for enhanced Atlantic Meridional Overturning Circulation (AMOC) compared to the present, but this hypothesis is not supported by the Pliocene Model Intercomparison Project (PlioMIP). Here, we adopt a new approach to assess variability in deep ocean circulation based on paired stable carbon (δ13C) and oxygen isotopes (δ18O) of benthic foraminifera, both (semi)conservative water mass tracers. Assuming that deep Atlantic benthic δ13C‐δ18O variability is predominantly driven by mixing, we extrapolate the δ13C‐δ18O data outside the sampled range to identify the end‐members. At least three end‐members are needed to explain the spatial δ13C‐δ18O variability in the deep North Atlantic Ocean: two Northern Component Water (NCW) and one Southern Component Water (SCW) water masses. We use a ternary mixing model to quantify the mixing proportions between SCW and NCW in the deep Atlantic Ocean. Our analysis includes new data from Ocean Drilling Program Sites 959 and 662 in the eastern equatorial Atlantic and suggests that the AMOC cell was deeper during the M2 glacial than during late Pliocene interglacials. Moreover, we identify a new cold and well‐ventilated water mass that was geographically restricted to the southeast Atlantic Ocean between 3.6 and 2.7 Ma and did not contribute significantly to the δ13C‐δ18O variability of the rest of the basin. This high‐δ13C high‐δ18O water mass has led to the misconception of a reduced latitudinal δ13C gradient. Our analyses show that the late Pliocene δ13C gradient between NCW and SCW was similar to the present‐day value of 1.1‰.


Introduction
The late Pliocene (3.6-2.6 Ma) interglacial climates may be the best past analogs for end 21st century climate (Haywood et al., 2011). Within this time frame, a particularly well-studied interval is the late Pliocene Mid-Piacenzian Warm Period (MPWP, 3.264-3.025 Ma). Global MPWP surface temperatures were on average 2-3°C higher than in the preindustrial, and a large part of this warming was concentrated at higher latitudes (Dowsett et al., 2010). Temperature proxies record peak Pliocene land surface warming of~19°C in the high Arctic (Ballantyne et al., 2010) and ocean surface warming of~8°C in the North Atlantic Ocean (Dowsett et al., 2016). Several mechanisms have been proposed to explain excessive warmth at high latitudes, including elevated atmospheric CO 2 levels (Lunt et al., 2012;Raymo et al., 1996), and strengthening of the Atlantic Meridional Overturning Circulation (AMOC) (Raymo et al., 1996;Zhang, Nisancioglu, & Ninnemann, 2013).
In previous studies, AMOC strength was typically assessed through spatial patterns in the stable carbon isotopic composition (δ 13 C) of deep sea benthic foraminifera, which is a proxy for the δ 13 C of dissolved inorganic carbon (δ 13 C DIC ). Within the major part of the modern Atlantic basin, δ 13 C DIC is not significantly impacted by organic matter remineralization, so that it generally acts as a conservative water mass tracer (Oppo & Fairbanks, 1987). This means that the δ 13 C DIC of a water mass will only change when it mixes with a water mass that has a different δ 13 C DIC signature. Consequently, deep Atlantic stable carbon isotope records reflect the site-specific relative contribution of 13 C-enriched deep water sourced from the North Atlantic and Arctic region (e.g., Northern Component Water [NCW]) to 13 C-depleted deep water sourced from the Southern Ocean and Antarctic regions (e.g., Southern Component Water [SCW]) (Kroopnick, 1980(Kroopnick, , 1985. Possible complications with using δ 13 C as a water mass tracer in palaeoceanographic reconstructions may arise when it loses its conservative properties due to low extremely low deep ocean ventilation rates (i.e., water mass aging) or local effects of remineralization in areas of high export productivity, for example, the Mackensen effect in the Southern Ocean (Mackensen et al., 1993).
Several authors have reported reduced spatial δ 13 C gradients in the Atlantic Ocean, which has been interpreted as evidence for increased NCW production and, consequently, for increased late Pliocene AMOC strength compared to the present (Haug & Tiedemann, 1998;Hodell & Venz, 2006;Poore et al., 2006;Ravelo & Andreasen, 2000;Raymo et al., 1992;Venz & Hodell, 2002). However, modeling experiments carried out in the context of the Pliocene Model Intercomparison Project (PlioMIP) fail to simulate consistently such an enhancement (Zhang, Nisancioglu, & Ninnemann, 2013). This apparent model-data discrepancy provokes further reflection on the interpretation of benthic isotope records.
The contribution of NCW and SCW to the deep Atlantic Ocean is commonly quantified using the binary mixing model that was formulated by Oppo and Fairbanks (1987) (Equations 1 and 2). We will refer to this model as the Binary Mixing Model (BMM).
This BMM is useful for determining a simple system but has potential limitations. First, the model requires knowledge of the isotopic composition of the northern-and southern-sourced deep-water end-members, which is problematic because the isotopic end-member may not always be perfectly preserved due to diagenesis of benthic foraminifera or not sampled by ocean drilling campaigns. Second, the model assumes that there are merely two end-members involved in deep ocean circulation, but like at present, there were likely multiple in the past (Bell, Jung, Kroon, Hodell, et al., 2015;Duplessy et al., 2002;Raymo et al., 2004).
For every benthic δ 13 C record there is an oxygen isotope (δ 18 O) counterpart, because these two variables are simultaneously measured in conventional carbonate stable isotope mass spectrometry. Foraminiferal δ 18 O reflects seawater δ 18 O (δ 18 O sw ) and temperature, which are both conservative water mass properties and can therefore be used to trace deep water masses (Bell et al., 2014;Bell, Jung, Kroon, Hodell, et al., 2015;Du et al., 2016;Duplessy et al., 2002;Raymo et al., 2004). Hence, complementary δ 18 O records may be extremely valuable in detecting different water mass end-members and quantifying past variability in deep ocean circulation, as shown by . We therefore propose to expand the BMM to account for both δ 13 C and δ 18 O variability in a ternary mixing model (TMM), similar to the work of Du et al. (2016) for the early Pleistocene.
The aim of our work is twofold. First, we present new late Pliocene δ 18 O and δ 13 C records from the deep eastern equatorial Atlantic. This largely neglected region is important to trace deep ocean circulation, because it connects the North Atlantic to the southeast Atlantic and is located upstream on the path of NCW to the Walvis Ridge. The past source of deep water to this region is not well understood Farmer et al., 2019). We generated a high-resolution stable isotope record for Ocean Drilling Program (ODP) Site 959 and tuned it to the LR04 benthic oxygen isotope stack (Lisiecki & Raymo, 2005) for accurate age control. We also present the late Pliocene δ 13 C record of the nearby ODP Site 662, located 1,100 km to the southwest from Site 959 ( Figure 1). Secondly, we reevaluate 15 previously published benthic isotope records together with our new records to identify the isotopic composition (paired the δ 13 C-δ 18 O) of the end-members that have filled the late Pliocene deep Atlantic (3.5-2.8 Ma, divided in seven time slices). We then use the TMM to quantify the mixing proportions of these end-members for a set of late Pliocene benthic isotope records from the north and tropical Atlantic. To assess if the results can be used for comparison to the modern, we compare the results of binary and ternary mixing, and compare Pliocene to Last Glacial Maximum (LGM), Holocene, and core top data.

Atlantic Ocean Deep Water Masses
The deep Atlantic is bathed by a mixture of NCW and SCW, which in the modern ocean are North Atlantic Deep Water (NADW) and Antarctic Bottom Water (AABW), respectively. These are, in turn, mixtures of several distinct water masses. NADW can be subdivided into a middle and deep component, North East Atlantic 10.1029/2019PA003804

Paleoceanography and Paleoclimatology
van der WEIJST ET AL.
The sourcing of Atlantic Ocean deep waters varies over time, in part because regions of deep-water formation may shift geographically in relation to, for example, sea ice conditions and surface ocean salinity and temperature, but also because the rates of deep-water formation may vary within each region. This variability occurs both on short (interdecadal) timescales (Ferreira & Kerr, 2017) and on geological timescales.
The spatial distribution of NCW and SCW in the deep Atlantic Ocean depends not only on water mass formation rates and densities, but also on basin bathymetry (Figure 1). The Mid-Atlantic Ridge (MAR) separates the eastern and western Atlantic basin at an average depth of 2,500 m. As a result, zonal transport of deep waters is largely limited to the larger fracture zones that punctuate the MAR, such as the Charlie-Gibbs Fracture Zone (CGFZ,~52°N, predominantly east-west flow), the Vema Fracture Zone (VFZ,~10°N, west-east flow), the Romanche Fracture Zone (RFZ, between~1°N and~1°S, west-east flow), and the Chain Fracture Zone (CFZ,~1°S, west-east flow). Another key feature is the northeast-southwest oriented Walvis Ridge, a topographic barrier that separates the Angola Basin from the Cape Basin in the southeast Atlantic. Restricted south-north flow occurs between 2,500 and 3,500 m, and abyssal flow is limited to two small passages that extent below 4,000 m, including the Walvis Passage (Connary & Ewing, 1974). The MAR and Walvis Ridge together impose a strong east-west differentiation of the South Atlantic Ocean, where the deep western basin is dominated by SCW and the eastern basin by NCW, although with 20-30% contribution of SCW contribution sourced through the RFZ and CFZ (Arhan et al., 2003;Bell et al., 2014).

Spatial δ 18 O and δ 13 C Variability
The chemical and physical properties of the NCW and SCW end-members, including salinity, temperature, δ 13 C DIC and δ 18 O sw , are largely determined by surface ocean conditions in regions of deep-water formation. The end-member δ 13 C DIC signal in source regions is predominantly a product of surface productivity and respiration, and sea water advection (Hodell & Venz, 2006;Ravelo & Hillaire-Marcel, 2007). The characteristic high-δ 13 C signature of NCW results from preferential fixation of the 12 C isotope during photosynthesis in the North Atlantic. Water mass aging is of minor influence on δ 13 C DIC in the Atlantic basin where residence times of water are relatively low compared to remineralization rates. Significant aging occurs when NCW is recirculated with Pacific and Indian Ocean deep water around the Antarctic continent in the Antarctic Circumpolar Current (ACC), giving rise to the characteristic low-δ 13 C DIC signature of the SCW that is formed around Antarctica. The magnitude of the aging effect in the Southern Ocean depends on ocean ventilation rates, which have been shown to vary as a function of the Antarctic front position and associated easterlies strength (Zhang, Nisancioglu, Chandler, et al., 2013). In the warm late Pliocene Southern Ocean, the polar front was likely at a more southerly position compared to the preindustrial, which intensified the easterly winds and decreased ocean stratification (Barron, 1996;McKay et al., 2012;Zhang, Nisancioglu, Chandler, et al., 2013). This may have reduced the NCW-SCW δ 13 C gradient compared to present.
Deep ocean foraminiferal δ 18 O is controlled by seawater temperature (δ 18 O decreases by~0.25‰/°C with increasing temperature (Shackleton, 1974)), and by the δ 18 O sw in source regions. Variability in δ 18 O sw reflects the balance between precipitation and evaporation, freshwater runoff, sea ice growth and melt processes and sea water advection, on top of global ice volume and the δ 18 O of the ice sheets on glacial-interglacial timescales (Frew et al., 2000;Ravelo & Hillaire-Marcel, 2007). Together, these processes determine the slope of the regionally and temporally variable δ 18 O sw -salinity (δ 18 O-S) relationship (Rohling & Bigg, 1998), which is traceable in the deep ocean. For example, in the modern ocean, NEADW and NWABW have distinct δ 18 O-S slopes before they mix to the relatively homogenous NADW (Frew et al., 2000). The slope of NEADW is less steep than that of NWABW as a result of the contribution of LSW to NWABW, which is characterized by a high salinity and relatively low δ 18 O due to regional sea ice formation (Frew et al., 2000). Although NEADW has a δ 18 O sw of +0.1‰ compared to NWABW, it is also 1-2°C warmer, which should result in an opposite isotopic gradient in calcite precipitated in these end-member water masses (Bell, Jung, Kroon, Hodell, et al., 2015;Frew et al., 2000). initially identified based on nannofossil biostratigraphy (Shin, 1998) and comprise nannofossil/foraminifer ooze. Recently, Vallé et al. (2017) generated a spliced composite record of the upper~100 m of Holes A, B, and C based on X-Ray Fluorescence scan data. They tuned Fe intensity from their new composite record to orbitally induced variations in insolation (sum of eccentricity, precession and tilt), derived from the La2004 astronomical solution (Laskar et al., 2004), constraining the late Pliocene interval with four cyclostratigraphic tie points. We established a revised age model of the upper Pliocene at Site 959 by visually correlating our new benthic δ 18 O record to the LR04 benthic isotope stack (Lisiecki & Raymo, 2005).

Site Descriptions and Age Models
ODP Site 662 (1°23.41′S, 11°44.35′W) was drilled during Leg 108 at a water depth of 3,820 m (Ruddiman et al., 1988) and is located at 1,100 km to the southwest of ODP Site 959 (Figure 1). The lithology consists mainly of nannofossil or nannofossil/foraminifer ooze. The age model of Site 662 was generated by Lisiecki and Raymo (2005) who used the δ 18 O record to construct the LR04 isotope stack.

Stable Isotope Analyses
For ODP Site 959, stable isotope analysis was performed on the epibenthic species Cibicides wuellerstorfi. Typically, three to six specimens were picked from the >212 μm size fraction. Calcite preservation was visually judged to be good to excellent throughout the interval, with many foraminifera being of a "glassy" appearance. Specimens were gently crushed, ultrasonically cleaned in ethanol to remove clays and coccoliths and subsequently rinsed with distilled water. Analyses were performed at two different setups at the Faculty of Geosciences, Utrecht. A Thermo-Finnigan Kiel III automated preparation system coupled to a Thermo-Finnigan MAT 253 mass spectrometer was used for 80 samples and a Thermo Finnigan GasBench-II carbonate preparation device coupled to a Thermo Finnigan Delta-V mass spectrometer was used for the remaining 162 samples. Due to the limited abundances of C. wuellerstorfi in our samples, we did not perform a cross-check on our own material. At the time of our measurements, both setups were running stably, and standards and replicate measurements performed on both systems showed consistent results. We find no systematic offsets between the results from the two instruments (supporting information Figure S2). For calibration to the international standard (Vienna Peedee belemnite [VPDB]), an in-house standard (Naxos for the Kiel and GasBench systems) and international standards were used (NBS-19 for the Kiel system and IAEA-CO-1 for the GasBench system). The analytical precision (standard deviation of the Naxos standard) of the Kiel system was better than 0.08‰ for δ 18 O and 0.04‰ for δ 13 C. On the GasBench system, analytical precision was better than 0.06‰ for δ 18 O and 0.04‰ for δ 13 C.
The oxygen and carbon isotope data for Site 662 were generated in the early 1990s at Massachusetts Institute of Technology using a VG Prism mass spectrometer with an analytical precision better than ±0.08‰ for δ 18 O and ±0.05‰ for δ 13 C. At this time, sample replicates from were better than ±0.3‰ (1σ) for both δ 18 O and δ 13 C (this includes inherent geologic variability). Foraminifera were picked from the >150 μm size fraction. Whenever possible, C. wuellerstorfi were analyzed; however, low abundances sometimes necessitated the analysis of Cibicides kullenbergi, Cibicides bradyii, or a mixture of Cibicides species. The use of species did not systematically impact the generated records ( Figure S3).
At both sites, the Cibicides δ 18 O data were adjusted with +0.64‰ to correct for vital effects between species and ambient seawater chemistry following Shackleton et al. (1995).

Treatment of Paired δ 13 C-δ 18 O Data
We follow a similar approach to Bell, Jung, Kroon, Hodell, et al. (2015), who visualize large quantities of data in comprehensive plots of paired δ 13 C-δ 18 O data, averaged over several Pliocene and Pleistocene time slices. We combine our new data with previously published high-resolution deep ocean stable isotope records from different ocean basins ( Figure 1). Age models that were originally constructed by tuning geochemical data to orbital solutions, or had asynchronous peak glacial δ 18 O values compared to the LR04 stack, were revised by visually aligning the δ 18 O data to the LR04 stack (i.e., Sites 925, 929, 659, 1148, and 1123, see also Figure S4). We initially defined seven separate time slices centered around the MPWP: three peak glacial stages (i.e., Marine Isotope Stage [MIS] M2, KM2, and G20), two phases of the MPWP (i.e., MPWP-1 and MPWP-2) and two bordering intervals ( Figure S4 and S5). For each interval, the δ 18 O and δ 13 C data were averaged. To assure data quality, averages estimated from fewer than five measurements or with very large standard errors were largely neglected in further analysis.
The majority of the data in our compilation were generated using tests of Cibicides wuellerstorfi, but in some records, other species were used as well. For details, the reader is referred to the original papers that are cited in the caption of Figure 1. To overcome species-specific offsets, all records were corrected following Shackleton et al. (1995).

Using Paired δ 13 C-δ 18 O Data to Identify Deep-Water End-Members
At present, calcite precipitated in the young NCW is characterized by relatively high δ 13 C and low δ 18 O, whereas calcite precipitated in the old and cold AABW has relatively low δ 13 C and high δ 18 O. Assuming that both δ 13 C and δ 18 O behave as conservative tracers within the Atlantic basin, the mixing of these two end-members can be expressed by a linear BMM in which the relative contribution of each water mass determines the δ 13 C-δ 18 O of foraminifera ( Figure 2a). However, the large spread in the Pliocene Atlantic Ocean δ 13 C-δ 18 O data indicates that a binary model is not sufficient to describe the mixing of deep Atlantic water masses (e.g., Bell et al., 2014;Raymo et al., 2004). Rather, a mixing scheme with not one but two NCW endmembers, one with relatively high δ 13 C-high δ 18 O and one with relatively high δ 13 C-low δ 18 O (Figure 2b), in addition to a low-δ 13 C high-δ 18 O SCW, could account for much of the spread seen in the north and tropical Atlantic sites.
SCW can be difficult to trace because enrichment in respired CO 2 leads to poor carbonate preservation. Previous studies used Pacific Ocean (e.g., Lang et al., 2016;Poore et al., 2006) or southeast Atlantic (e.g., Venz & Hodell, 2002) δ 13 C records to represent the SCW end-member. However, the Pacific Ocean records are probably substantially affected by aging, which lowers their δ 13 C signature. This may be corrected for using the modern offset between Southern Ocean and Pacific Ocean deep waters of~0.5‰ (Kroopnick, 1985), but this approach is problematic with changing modes of ocean circulation through time. Moreover, a recent study demonstrated that during the late Pliocene, deep waters may have formed in the Subarctic Pacific Ocean (Burls et al., 2017). This Pacific Meridional Overturning Circulation (PMOC) cell may have further obscured the original δ 13 C-δ 18 O signature of SCW in the equatorial and north Pacific.
In contrast to many previous studies where a specific isotopic record was assigned an end-member status, we choose to estimate the end-member values from the spread of δ 13 C-δ 18 O data from the north and tropical Atlantic, Indian Ocean Site 758 and south Pacific Site 1123, assuming that the majority of the spread in this data is explained by NCW and SCW mixing. Using the MinBoundSuite algorithm of D'Errico (2020), we calculate a range of possible ternary mixing envelopes for each time slice, accounting for the error of the δ 13 C and δ 18 O estimates; see also Figure S7 and Text S3.

Assessing Water Mass Mixing With the TMM
The end-member δ 13 C-δ 18 O estimates that result from our analyses are used in both the BMM (Equation 1) and TMM to quantify NCW and SCW in the deep tropical and North Atlantic. The TMM describes the δ 13 C-δ 18 O data at a particular site as a function of ternary mixing between one SCW and two NCW end-members (Figure 2b). Mixing is described in two steps: First, NCW1 mixes with NCW2 to form a site-specific NCW′ and then NCW′ mixes with SCW. This two-step mixing is preferred over equal mixing between the three end-members because we assume that mixing between the two NCW end-members precedes mixing between NCW and SCW due to the geographic placement of the source regions, analogous to modern NADW. The relative contribution of NCW1 and NCW2 to NCW′ is described by a linear mixing function, as is the relative contribution of SCW and NCW′ to a specific site in the deep Atlantic Ocean. The equations of %SCW, %NCW1, and %NCW2 (Equation 3-5, respectively) require the calculation of four dimensionless parameters: a is the length between NCW1 and NCW′, b is the length between NCW2 and NCW′, c is the length between the sample and NCW′, and d is the length between the sample and SCW. All calculations are described in the supporting information (Text S1).
To account for the uncertainty in the size and shape of the mixing envelopes and the exact position of each site in δ 13 C-δ 18 O space, we performed 1,000 Monte Carlo iterations on the Pliocene data set to calculate water mass mixing proportions, see also Text S3.
To our knowledge, there is no comprehensive Holocene paired deep ocean benthic foraminifer δ 13 C-δ 18 O compilation available to evaluate our method for the present day. We therefore asses if the compilation of core top data of Schmittner et al. (2017) might serve as a suitable alternative. Application of the TMM and BMM to core top data ( Figure S1) produces nearly identical mixing proportions of NCW and SCW, although the relative volumes of SCW are slightly higher at high %SCW sites according to the TMM model ( Figure 3). However, both models produce much higher %SCW estimates compared to models that depend on physical and chemical tracers measured in deep ocean water rather than sediment (Ferreira & Kerr, 2017;Johnson, 2008). We further test the TMM with the LGM data set of (Duplessy et al. (2002) in Text S4 of the supporting information.

Correlation of Site 959 to the LR04 Stack
Alignment of our new ODP Site 959 δ 18 O data with the LR04 chronology ( Figure 4b) shows a high degree of correspondence, both regarding absolute values as well as the amplitude and timescale of variability, and hence provides a robust framework for correlation to other deep ocean sites. The initial tuning of X-ray fluorescence (XRF) data to the Eccentricity-Tilt-Precession (ETP) curve (Vallé et al., 2017) corresponds also very well but inconsistencies in the initial tuning and/or the adopted phase relationships between the ETP target curve and the XRF data led to errors of up to 50 kyr for some intervals (Figure 4a). The investigated depth interval at Site 959 spans an age of~2.80 to 3.50 Ma. Sedimentation rates, calculated from the slightly stretched and compressed rmcd depth scale, vary between 1.0 and 2.8 cm/kyr around a mean of 1.8 cm/ kyr (Figure 4c). The good correspondence with the LR04 record suggests continuous sedimentation and a negligible effect of core gaps or duplicated sections.

Variability in Late Pliocene Eastern Equatorial Atlantic Stable Isotope Records
We compare our new isotope records of ODP Site 959 with the published δ 18 O (Lisiecki & Raymo, 2005) Table S1 for age bounds.

Paleoceanography and Paleoclimatology
van der WEIJST ET AL.
therefore is expected to be 18 O-enriched due to cooler and/or more saline waters. We cannot explain this reversed gradient with diagenesis or interlaboratory offsets (Text S2 and Figure S6). Figure 6 shows the distribution of the paired isotope data in δ 13 C-δ 18 O space. The estimated end-member δ 13 C-δ 18 O ranges for various late Pliocene time intervals are listed in Table 1. The δ 13 C-δ 18 O data are largely clustered by geography, separating the various ocean basins. These patterns are similar to those seen in the modern ocean, consistent with previously work (Bell, Jung, Kroon, Hodell, et al., 2015). Our compilation shows that the southeast Atlantic record from Site 704 is inadequate to represent the SCW end-member because its δ 13 C signature is higher than that of several north and tropical Atlantic sites ( Figure 6). Southeast Atlantic Site 704 requires an additional end-member to explain their remarkable high δ 13 C-high δ 18 O signature, as do Walvis Ridge Sites 1264 and 1267 in some time slices (Figure 2c). This will further be discussed in section 5.2.

Water Mass Mixing in the Late Pliocene Deep Atlantic Ocean
Instead, we assume that Indian Ocean Site 758 and south Pacific Site 1123 best approximate the SCW δ 13 C-δ 18 O. Both sites record a high δ 18 O-low δ 13 C signature, which is in line with present-day AABW

Paleoceanography and Paleoclimatology
van der WEIJST ET AL. (Bell, Jung, Kroon, Hodell, et al., 2015). In analogy with the present, the δ 13 C gradient between these sites and the tropical Pacific indicates a minor aging effect on these records. Also, both sites would have been unaffected by the potential presence of Pliocene PMOC (Burls et al., 2017).
To explain all δ 13 C-δ 18 O variability in the tropical and North Atlantic sites with a ternary mixing scheme, the mixing plane has to be extrapolated beyond currently available isotopic data. The convex hull outlines represent mixing lines between two of the end-members and data that plots on the plane represents mixing of three end-members ( Figure 6). We recognize two distinct NCW end-members, NCW1 with a high-δ 13 C high-δ 18 O signature, and NCW2 with a high-δ 13 C low-δ 18 O signature and a SCW end-member with a low δ 13 C high-δ 18 O signature (Table 1).
In the G20 time slice, the convex hull is largely determined by the δ 13 C-δ 18 O of Site 925, which does not meet the data quality criteria defined in section 3.3. Similarly, the shape of the KM2 mixing plane does not resemble the other time slices, possibly because of undersampling of the spatial δ 13 C-δ 18 O variability during this time slice. Because these issues impose large uncertainty on the following reconstructions, mixing patterns were not further assessed for the G20 and KM20 time slice.
Using currently available data in the TMM, we were able to make a detailed reconstruction of water mass mixing proportions during the M2 glacial stage, an early and late MPWP time slice, and of the background state during two bordering intervals (Figures 7 and 8 and Table S2). The uncertainties in the reconstructed water mass proportions are largest for %NCW1 and %NCW2 (see Table S2), which means that it is difficult to completely deconvolve these signals, particularly during the M2 glacial. The larger uncertainty in this particular time slice is predominantly a result of the smaller time window, and consequently a lower number of samples used in the compilation. Moreover, the NCW1 and NCW2 end-members are relatively closely spaced during M2, which makes it more difficult to deconvolve the signals than in the other time slices.
The uncertainties of %SCW are relatively low, which suggests that the model performs well in deconvolving northern-and southern-sourced water masses.
Figures 7 and 8 reveal striking temporal variability of water mass mixing during the late Pliocene. Most sites record a mixed signal of NCW1, NCW2, and SCW. None of the sites record a clear linear trend in the relative contribution of any of the water masses, indicating that processes acting on shorter timescales were more important in shaping late Pliocene thermohaline circulation than processes acting on longer timescales,  Table S2 for numerical results). Relative abundances were also calculated for the South Atlantic sites when their data would fit within the mixing envelope (low-opacity pie charts). Because they are influenced by an additional water mass (SEADW, see section 5.2), the TMM may be inappropriate to deconvolve the isotope records from these sites. The bathymetric profile was generated in Ocean Data View (Schlitzer, 2019).

Paleoceanography and Paleoclimatology
van der WEIJST ET AL.

Source Regions and Circulation Patterns
Although the absolute Pliocene %SCW estimates may be too high (see section 3.5), the spatial distribution of SCW can very well be explained by modern circulation patterns. SCW generally shows an increasing abundance with increasing depth and increasing proximity to the southern hemisphere (Figures 7 and 8 for the Pliocene and Figures S8 and S9 for the LGM). At present, the majority of SCW (~1.1 Sv) enters the eastern North Atlantic via fracture Zones in the MAR and~0.3 Sv SCW ends up in the western North Atlantic basin (Rhein et al., 1998). Similarly, during the late Pliocene, SCW was more abundant in the eastern than the western Atlantic basin, especially at shallower depths. Notably, the interface between SCW and NCW was much sharper in the western than in the eastern basin. In the western basin, the most northerly site (Site U1313 at 41°N, 3,430 m) records a smaller influence of SCW compared to the most northerly site in the eastern basin (Site 610 at 53°N, 2,420 m).
The high-δ 13 C intermediate-δ 18 O NCW1 was typically more abundant at shallower depths than the high-δ 13 C low-δ 18 O NCW2. During M2, stratification was at a maximum in the Atlantic basin and NCW1 and NCW2 did not mix as thoroughly as in the other time slices (see also section 5.1.2). At this time, NCW2 was highly abundant in the eastern, but not in the western basin, a pattern that is best explained by an upstream source of NCW2 in the eastern basin. Together, these observations suggest that NCW1 was a Pliocene equivalent of the modern LSW, sourced from the northwest Atlantic, and NCW2 of the denser overflow waters from the Nordic Seas (Figure 1).
While NCW1 was more abundant at shallower depths, that is, less dense than NCW2, although the reversed vertical δ 18 O gradient, which is well expressed in the eastern equatorial Atlantic, suggests otherwise. Assuming the present-day~1°C temperature gradient between Sites 662 and 959 (Locarnini et al., 2013) and the modern NADW δ 18 O-S slope of 0.61Δδ 18 O/salinity unit (Frew et al., 2000), NCW1 must have had a 0.7-0.9 higher salinity than NCW2. However, this would likely have resulted in a physically unattainable vertical density gradient. Instead, the reversed δ 18 O gradient is more likely the result of a different δ 18 O-S relationship in the shallower NCW1 and deeper NCW2: the deeper NCW2 must have had a steeper δ 18 O-S slope than NCW1. Next, we will briefly review the possible flow patterns for each time slice.

From 3500 ka to M2
In the 3500 ka to M2 time slice (Figures 7 and 8, first panels), the relative water mass distribution in the North Atlantic can be explained by a similar-to-present-day circulation pattern, with the largest exception being that SCW reached higher northern latitudes and was not limited to abyssal depths, but rather mixed with NCW1 and NCW2 at shallower sites such as Site 959 (2,090 m). At Site 925 (3,040 m), NCW1 and NCW2 mixed in roughly equal proportions. However, NCW1 was absent at Site 929 (4,360 m), suggesting substantial density stratification within the Deep Western Boundary Current (DWBC). Downstream from Site 925, at Sites 662 (3,820 m) and Site 659 (3,080 m) across the RFZ and CFZ, %NCW2 was much higher than % NCW1. This suggests the presence of an additional pathway for NCW2, on the east side of the MAR. This branch must have been deeper than Site 959, which does not seem affected by this current. The δ 13 C-δ 18 O signal at Sites 1264, 1267, and 704 is best explained by the presence of a fourth water mass, which we will refer to as South East Atlantic Deep Water (SEADW) after its locality. It may have been the Pliocene analog to Antarctic Intermediate Water and will further be discussed in section 5.2.

M2 Glacial
During the M2 glacial (Figures 7 and 8, second panels), mixing between NCW1, NCW2 and SCW was at a minimum, resulting in a relatively stratified Atlantic Ocean. SCW filled large parts of the abyssal Atlantic, as recorded at Site 929. Slightly shallower, at Site 662, there was approximately an equal mix between SCW and NCW2, the densest of the NCW masses. As in the 3500 ka to M2 time slice, there seems to have been a NCW2 current flowing along the eastern flank of the MAR, explaining the high %NCW2 recorded at Site U1308 (3,900 m). The absence of NCW2 at Site 610 (2,420 m) indicates that the main body of this water mass resided deeper. NCW1 generally dominated the shallowest sites, and at intermediate depths, NCW2 was more abundant in the eastern basin and NCW1 in the western basin, as indicated by the relative water mass abundances at Sites 925 and 659. It remains unclear if the ternary mixing scheme we used for the tropical and North Atlantic sites is also appropriate for Site 1264, but if so, then there must have been a substantial supply of NCW1 through the RFZ and CFZ and the export of SEADW must have been reduced.

MPWP-1
The water mass distribution during the MPWP-1 time slice (Figures 7 and 8, third panels) was similar to the 3500 ka to M2 time slice, with a few notable exceptions. First of all, SCW was less abundant in the Atlantic, especially at shallower depths (e.g., Site 959) and high northern latitudes (e.g., Site U1313). Moreover, abyssal Site 929 records a high abundance of NCW1, which is unexpected because of its low density. Furthermore, this NCW1 distribution is difficult to reconcile with the other sites in the Atlantic. A possible explanation to this aberrant profile is that SEADW was more prominent in this time slice and reached the deeper Ceara Rise Site 925, thereby lowering the δ 18 O signature toward NCW1 values. Just as in the 3500 ka to M2 time slice, Site 704 and both Walvis Ridge sites fall well outside the ternary mixing scheme, suggesting a strong, albeit local, influence of SEADW. At Site 925, NCW1 makes up 40% of the total NCW. Downstream, at Site 662 across the RFZ and CFZ, the ratio between NCW1 and NCW2 is identical, suggesting that the NCW2 branch on the eastern side of the MAR did not reach as far south as during the previously discussed time slices. Still, its presence during this time slice is indicated by the high %NCW2 at Site 662.

MPWP-2
In the MPWP-2 time slice (Figures 7 and 8, fourth panels), the ratio between NCW1 and NCW2 is similar at Sites 925 and 662, but at Site 659, NCW1 is much more abundant. This suggests that the NCW2 branch on the east side of the MAR weakened or ceased to exist. This is further corroborated by the high abundance of NCW1 instead of NCW2 at Site U1308. One explanation for these patterns is reversed flow through the CGFZ at 52°N, potentially forced by North Atlantic Current (NAC) interference with deep water. At present, the eastward flowing NAC crosses the MAR between 45°N and 53°N, with the main branch typically positioned just south of the CGFZ (Bower & Von Appen, 2008). A temporary northward excursion of the NAC can reverse the flow direction in the CGFZ, effectively blocking the westward flowing ISOW (Bower & Von Appen, 2008;Schott et al., 1999). At present, this occurs at 0-35% of the time (Bower & Von Appen, 2008). Could a prolonged northward shift of the NAC in the Pliocene raise this number?
Prior to M2, the NAC was intensified and its effects were traceable at Site 982 (57°N), but several lines of evidence suggest it weakened and shifted southward over the course of the late Pliocene (Karas et al., 2020; 10.1029/2019PA003804

Paleoceanography and Paleoclimatology
van der WEIJST ET AL. Lawrence et al., 2009;Naafs et al., 2010), reaching Site U1313 at 41°N around~2.6 Ma (Bolton et al., 2018;Hennissen et al., 2014). We speculate that the NAC was positioned over the CGFZ during MPWP-2, blocking westward NCW2 flow and promoting eastward NCW1 flow. If the westward flow of NCW2 via the CGFZ was largely or completely blocked, the westward flow in the Denmark strait may have intensified to compensate. Another striking feature during MPWP-2 is the high abundance of SCW at Site 610, which may point to an alternative pathway for northward SCW flow, such as a crossing of the MAR at the VFZ. However, more data is needed to constrain circulation patterns in this much detail. Alternatively, although we have no concrete indications for it, this exceptional reconstruction may be affected by calcite preservation or temporary changes in seawater chemistry.
The δ 13 C-δ 18 O signature at both Walvis Ridge Sites falls within the ternary mixing scheme, which theoretically enables the deconvolution of the signal into NCW1, NCW2, and SCW contribution. The reconstructed mixing proportions (Figure 7) at Site 1264 (2,510 m) are similar to the M2 time slice and seem realistic if these were indeed the water masses that mixed at this site. However, δ 13 C-δ 18 O signature is identical at these sites, which leads to an identical water mass reconstruction at Site 1267 (4,360 m): the low-density NCW1 would have been dominant and the dense SCW would have been present in modest quantities.
Considering the depth at Site 1267 and its distance to the NCW source regions, this reconstruction seems unrealistic. The δ 13 C-δ 18 O signature at both sites may rather be shaped predominantly by SEADW, with a contribution of NCW1 and/or NCW2 ( Figure 6, orange convex hull). In contrast to the MPWP-1 time slice, there is no evidence of SEADW at Site 929.

G20 to 2800 ka
In the G20 to 2800 ka time slice (Figures 7 and 8, fifth panels), the mixing proportions are roughly similar to the MPWP-2 time slice. The biggest change relative to the previous time slice occurred at Site 610, which returns to a state with mixing proportions most similar to the 3500 ka to M2 time slice, with a dominant NCW1 component mixed with NCW2. Like in the MPWP-2 time slice, Site U1308 seems to record a signal that is resembles the western Atlantic basin, suggesting persistent eastward flow of NCW1 through the CGFZ, potentially influenced an intensified NAC at a higher latitude. Walvis Ridge Site 1264 may be recording an increased NCW2 influence compared to earlier time slices, although it remains unclear if the data are suited to be used in our reconstructions, because they may be biased by SEADW.

SEADW and the Atlantic δ 13 C Gradient
The majority of the late Pliocene Atlantic δ 13 C-δ 18 O data can be explained with the mixing of NCW1, NCW2 and SCW. However, Site 704 has an unusual high-δ 13 C high-δ 18 O signal compared to the other sites ( Figure 6), as do Walvis Ridges Sites 1264 and 1267 in certain time slices. Therefore, the mixing scheme requires an additional end-member, which we name SEADW. The hypothesized mixing scheme for Sites 704, 1264, and 1267 is illustrated in Figure 6. This interpretation differs substantially from previous work in which only the latitudinal δ 13 C gradients were considered and a high δ 13 C signature in the southeast Atlantic was explained by a strong contribution of NCW (Hodell & Venz, 2006;Patterson et al., 2018). However, the high δ 18 O values in the late Pliocene southeast Atlantic clearly show a different origin for these waters.
It seems that the spatial extent of SEADW was limited. The Walvis Ridge is an effective geographic barrier, although minor volumes of deep water cross it via the Walvis Passage (Connary & Ewing, 1974). This southnorth connection could explain the similarity between isotopic values at Site 1267 on the north face of the Walvis Ridge, Site 1264, located on top of the ridge, and Site 704, south of the ridge. Sites 662 and 959 show no signs of SEADW contribution, indicating that SEADW did not reach other sites beyond the Walvis Ridge, except maybe Site 929 during the MPWP-1 time slice (see section 5.1.2). SEADW may also be traceable in the Pacific basin, but data from the South Pacific is currently too scarce.
Compared to the SCW end-member, δ 13 C is~0.5‰ enriched at Site 704, and, except during the M2 glacial stage, δ 18 O is similar or higher. Together with the apparently limited geographic extent of SEADW, this isotopic signature suggests that it could be a well-ventilated southern-sourced water mass. It was perhaps similar to the modern Circumpolar Deep Water (CDW), which forms from vigorous mixing of both southern-and northern-sourced water masses in the ACC (Johnson, 2008;Ullermann et al., 2016). Indeed, also in the modern ocean, calcite precipitated from CDW is enriched in 13 C compared to calcite precipitated in 10.1029/2019PA003804

Paleoceanography and Paleoclimatology
van der WEIJST ET AL.
From a long-term perspective, the best diagnostic for the presence of the relatively cool and well-ventilated SEADW may be a positive δ 13 C gradient between Site 704 and SCW-dominated Site 1123, which occurs between 3.6 Ma and 2.7 Ma (Patterson et al., 2018). Before 3.6 Ma, both δ 13 C and δ 18 O were similar at Site 704 and Site 1123, which suggests no distinct SEADW. With reduced ice sheets and warmer Southern Ocean temperatures in the Early Pliocene (McKay et al., 2012;Patterson et al., 2018), the waters formed from the ACC may not have been dense enough to sink to the deep ocean. At 2.7 Ma, the δ 13 C gradient between Site 704 and Site 1123 reversed, possibly related to reduced ventilation of the Southern Ocean caused by increased seasonal sea ice (Hillenbrand & Cortese, 2006;McKay et al., 2012;Patterson et al., 2018), in combination with a southward shift of the polar front and a decrease in the easterlies strength (Zhang, Nisancioglu, Chandler, et al., 2013). From this point in time, SEADW can no longer be recognized as a high-δ 13 C high-δ 18 O southern-sourced water mass, although an analog water mass may still have formed by sinking waters in the ACC.
In previous work, the reduced δ 13 C gradient between the North Atlantic and the South Atlantic was ascribed to an strong influence of NCW on the South Atlantic, linked to an extensive and/or strong AMOC cell (Lang et al., 2016;Poore et al., 2006;Ravelo & Andreasen, 2000;Raymo et al., 1992;Venz & Hodell, 2002). However, our compilation shows that existing isotope records from the late Pliocene southeast Atlantic do not represent the SCW end-member but rather reflect a geographically limited intermediate to deep water mass (SEADW). Assuming that Indian Ocean and south Pacific δ 13 C values are a better approximation of the SCW δ 13 C, it seems that a pure SCW signal has not yet been sampled in the South Atlantic ( Figure 6). In fact, the average late Pliocene δ 13 C gradient of 1.05 (±0.10‰ SE) between our SCW estimate and the North Atlantic is identical to the present value of 1.1‰ in deep Atlantic core tops (Schmittner et al., 2017; Figure 3), although it decreased with >0.25‰ between 3.5 and 2.8 Ma ( Figure 6 and Table 1).
The reduced δ 13 C gradient between the North Atlantic and the southeast Atlantic has alternatively been interpreted in terms of increased ventilation (Hodell & Venz, 2006). During the late Pliocene, the polar front was at a more southerly position compared to the preindustrial, which would have intensified the easterly winds and decreased upper ocean stratification, leading to increased ventilation and a lower 13 C signature of southern-sourced deep waters (Barron, 1996;Zhang, Nisancioglu, Chandler, et al., 2013). Interestingly, simulations of the NorESM-L GCM described by Zhang, Nisancioglu, and Ninnemann (2013) show that intermediate depths (2-4 km) in the Southern Ocean were occupied by a geographically restricted well-ventilated cool water mass with a relatively low salinity. They also show that the denser, more saline waters below 4,000 m were poorly ventilated. This is in agreement with the contrasting isotopic signatures in the southeast Atlantic (high-δ 13 C high-δ 18 O) on the one hand and the South Pacific and Indian Ocean (low-δ 13 C high-δ 18 O) on the other hand. Our reconstructions show that this denser SCW, not the shallower, well-ventilated SEADW, was the bulk of southern-sourced water exported to higher northern latitudes. Hence, increased ventilation of the Southern Ocean did indeed occur but it did not drive a basin-scale reduction of the latitudinal δ 13 C gradient.

Binary Versus Ternary Mixing: Comparing Past and Present
The BMM and TMM produce nearly identical %SCW estimates from Atlantic core top data (Figure 3), which suggests that they perform equally well under the condition that the full deep ocean δ 13 C range is sampled. However, because this condition is typically not met in paleostudies, complementing δ 13 C with δ 18 O data may help to better identify deep-water end-members ( Figure 6). Moreover, the TMM allows for a more detailed reconstruction of the mixing of the highly variable NCW sources.
Because both the BMM and TMM overestimate %SCW in the North Atlantic when applied to core top data, caution is warranted when interpreting the reconstructed mixing volumes at face value. Although the Pliocene %SCW estimates seem high, they typically do not exceed the upper range of core top and Holocene δ 13 C-based estimates ( Figure S10), suggesting that the mixing volumes between NCW and SCW were similar to present day. Interestingly, a TMM-based reconstruction of %SCW in the LGM north Atlantic shows similar results ( Figure S10). This contradicts the traditional δ 13 C-based hypothesis that SCW dominated the LGM deep Atlantic (e.g., Boyle & Keigwin, 1982;Ferrari et al., 2014;Negre 10.1029/2019PA003804 et al., 2010. Instead, it supports the more nuanced view presented in recent papers (Du et al., 2020;Gebbie, 2014;Howe et al., 2016;Pöppelmeier et al., 2020), where low δ 13 C-values in the North Atlantic are explained by reduced ventilation and sluggish circulation. However, even if the additional δ 18 O tracer may improve water mass reconstructions, the TMM estimates are still partially determined by δ 13 C variability, and sluggish circulation might have interfered with the conservative properties of δ 13 C during the LGM.
Overturning rates during the Pliocene M2 glacial are poorly constrained, but because the LGM was a much stronger glacial than M2, with 3-10 times as much ice volume accumulation (Berends et al., 2019), we expect the magnitude of the aging effect to be much smaller.
In our reconstructions, northern-sourced waters dominated the deep Atlantic during the M2 glacial, except for the low latitude abyssal zones where southern-sourced water was abundant (Figures 7, 8, and S10). The AMOC cell was relatively deep and density stratified. With the transition into the warmer late Pliocene interglacial states, AMOC shallowed ( Figure 9). Interestingly, PlioMIP simulations show that high latitude warming can be sustained without significant changes in AMOC and ocean heat transport compared to the preindustrial (Zhang, Nisancioglu, & Ninnemann, 2013). However, considering the stark contrast in water mass distribution between M2 and the late Pliocene interglacials, we suggest that the relationship between AMOC and Pliocene climate variability deserves further attention.
Finally, we here present variations in deep Atlantic Ocean circulation in the Pliocene but further analyses are required for a proper comparison between the past and present. First of all, it is important to note that the core top data are biased heavily toward the east Atlantic ( Figure S1). It also lacks the time domain integrated in our Pliocene time intervals, which implies that it is influenced by variability on short (interdecadal) timescales (e.g., Ferreira & Kerr, 2017). Moreover, because in this work all data are equally weighed in the calculation of the mixing envelope, the end-member reconstructions are sensitive to outliers. In order to better compare the present and past situations, and to further validate the TMM, it is important to include more core top data from the west Atlantic and to critically define data selection criteria. For quantitative comparison of paleodata to the modern situation, we also recommend evaluating the TMM using a comprehensive late Holocene paired δ 13 C-δ 18 O compilation in future work.

Conclusions
We analyzed coupled δ 13 C-δ 18 O data in a compilation of 16 late Pliocene benthic isotope records, including new data from eastern equatorial Atlantic Sites 959 and 662. This compilation shows that the deep Atlantic Ocean was bathed by a mixture of at least four end-member water masses, of which at least three contributed significantly to the large spread in δ 13 C-δ 18 O in the deep north and tropical Atlantic. These comprise one SCW (low-δ 13 C high-δ 18 O) and two NCW end-members, of which one was formed in the Labrador Sea region (NCW1; high-δ 13 C intermediate-δ 18 O) and a denser one in the Nordic Seas (NCW2 high-δ 13 C low-δ 18 O). NCW1 and NCW2 had different δ 18 O-S relationships; the slope of NCW2 was relatively steep compared to NCW1, which caused a reversed vertical δ 18 O gradient in the eastern equatorial Atlantic. Figure 9. Simplified illustration of the water mass distribution in the Atlantic basin during the M2 glacial and a generalized state of the late Pliocene interglacials, which were characterized by a shallower AMOC. During interglacials, SCW was more abundant in the Atlantic basin; more mixing occurred between NCW1, NCW2, and SCW.
The fourth high-δ 13 C high-δ 18 O end-member, which we term SEADW, was geographically restricted to the southeast Atlantic Sites. Compared to NCW1 and NCW2, SEADW was a relatively cool and well-ventilated southern-sourced water mass, likely analogous to modern CDW. It was prominent in warm interglacial stages between 3.6 and 2.8 Ma. The high SEADW δ 13 C signature previously led to overestimates in late Pliocene latitudinal deep ocean δ 13 C gradients. Based on δ 13 C values of Indian and South Pacific Ocean sites, we conclude that the δ 13 C gradient in the late Pliocene interglacial deep Atlantic was similar to the present value of 1.1‰.
Based on a TMM that accounts for the spatial variability in both δ 13 C and δ 18 O, we show that the distribution of NCW1, NCW2, and SCW was different at both sides of the MAR. In the western basin, NCW1 overlay the denser NCW2 and the deeper parts bathed in SCW. Except for during the M2 glacial, the three water masses were better mixed in the eastern basin, possibly due to water transport from the western to the eastern basin through the Vema, Romanche and Chain fracture zones. During the M2 glacial, the Atlantic was relatively stratified compared to the Pliocene interglacials. SCW was highly abundant at the deepest tropical sites, but NCW1 and NCW2 dominated the deep North Atlantic. The Pliocene interglacial NCW and SCW mixing proportions fall largely within the modern spatial mixing variability as inferred from core top and Holocene data. However, the absence of a late Holocene paired δ 13 C-δ 18 O data set limits the possibility to investigate site-to-site differences between past and present.
Finally, our analyses of relatively short late Pliocene time slices reveal a vast variety of mixing patterns. This indicates that thermohaline circulation varied strongly on <100 kyr timescales, likely due to highly variable surface ocean conditions in source regions and interplay with surface currents such as the NAC. Our understanding of ocean circulation changes on shorter timescales may therefore benefit from better quantifications of the variability of deep-water mixing in the past and the TMM is a useful tool for this.