Pliocene decoupling of equatorial Pacific temperature and pH gradients

Ocean dynamics in the equatorial Pacific drive tropical climate patterns that affect marine and terrestrial ecosystems worldwide. How this region will respond to global warming has profound implications for global climate, economic stability and ecosystem health. As a result, numerous studies have investigated equatorial Pacific dynamics during the Pliocene (5.3–2.6 million years ago) and late Miocene (around 6 million years ago) as an analogue for the future behaviour of the region under global warming1–12. Palaeoceanographic records from this time present an apparent paradox with proxy evidence of a reduced east–west sea surface temperature gradient along the equatorial Pacific1,3,7,8—indicative of reduced wind-driven upwelling—conflicting with evidence of enhanced biological productivity in the east Pacific13–15 that typically results from stronger upwelling. Here we reconcile these observations by providing new evidence for a radically different-from-modern circulation regime in the early Pliocene/late Miocene16 that results in older, more acidic and more nutrient-rich water reaching the equatorial Pacific. These results provide a mechanism for enhanced productivity in the early Pliocene/late Miocene east Pacific even in the presence of weaker wind-driven upwelling. Our findings shed new light on equatorial Pacific dynamics and help to constrain the potential changes they will undergo in the near future, given that the Earth is expected to reach Pliocene-like levels of warming in the next century. New proxy data for ocean pH and an ocean–atmosphere model show that a radically different ocean circulation led to decoupling of ocean productivity and upwelling in the equatorial Pacific Ocean 3–6 million years ago.

A salient feature of the modern tropical Pacific climate is a pronounced zonal (east-west) sea surface temperature (SST) gradient along the equator, with mean temperature decreasing from around 29 °C in the western Pacific warm pool to around 23 °C in the eastern Pacific equatorial cold tongue 2 . This zonal SST gradient is tightly coupled to the Pacific zonal atmospheric circulation-known as the Walker cell-whose winds in turn control the zonal tilt of the oceanic thermocline and the strength of equatorial upwelling via the Bjerknes feedback 17 . The Pacific zonal SST gradient thus plays a key role in defining the climate of the tropics. Variations in this gradient produce the El Niño-Southern Oscillation-a dominant mode of interannual climate variability in the tropics and beyond 18 . In addition, vigorous upwelling in the eastern tropical Pacific brings low-pH, nutrient-rich water up into the surface ocean, supporting highly productive marine ecosystems that are important to the economic health and food security of adjacent nations 19,20 .
Despite the importance of the Pacific zonal SST gradient for climate and marine ecosystems, considerable uncertainty persists in how this ocean feature will respond to global warming 21 . One approach is to turn to past warm climates such as the Pliocene epoch (5.3-2.6 million years ago (Ma))-the last time in Earth's history when atmospheric CO 2 concentrations were similar (300-450 ppm) to the anthropogenically forced levels seen today 10,22 . By investigating Earth's climate during past warm periods such as the Pliocene (~2-4 °C warmer temperatures than modern [23][24][25][26][27][28] ) and indeed even the latest Miocene (~6 Ma, ~5 °C warmer 29 ), it is possible to develop a better understanding of tropical Pacific dynamics during warmer near-equilibrium climate states.
Existing records of the early Pliocene and late Miocene equatorial Pacific present a paradox. Today, a decrease in the zonal temperature gradient corresponds to reduced biological productivity in the eastern equatorial Pacific owing to reduced wind-driven upwelling of nutrient-rich waters-as occurs during an El Niño event. Perplexingly, records from the Pliocene and late Miocene seem to show a decoupling between the zonal SST gradient and biological productivity-with evidence for decreased zonal SST gradients 1,3,7,8 but increased productivity [13][14][15] . One solution to this paradox would be for the upwelled water to have been more nutrient-rich relative to modern, such that a smaller volume of upwelled water would support greater productivity.
Here we use the boronpH proxy to reconstruct the evolution of the Pacific zonal pH gradient along the equator since the late Miocene to test this idea that circulation changes delivered older, more nutrient-rich water to the eastern equatorial Pacific (EEP). This test relies on the fact that deep water that is upwelled to the surface contains https://doi.org/10.1038/s41586-021-03884-7 Article more nutrients than surface water owing to the remineralization of organic matter, which liberates nutrients along with CO 2 , resulting in a decrease in the upwelled water's pH. As a result, the upwelling of older water masses contributes proportionately more acidity and nutrients to the surface, owing to their prolonged accumulation of respired organic matter. There are four primary, but not necessarily mutually exclusive, ways to upwell relatively older water into the EEP: (1) the water mass source could be further away, (2) the path to the EEP could be longer, (3) the flow could be slower, and/or (4) the water mass could upwell from deeper depths with relatively older water.
We focus our efforts on pH as there is an ongoing debate as to the magnitude of the zonal temperature gradient reduction at the time, with estimates for the reduction ranging from around 1 °C (refs. 3,6,11 ) to around 3.5 °C (refs. 1,4,10,12,30,31 ) relative to modern, depending on the proxy, calibration, sites used and time interval of focus. Although the magnitude is debated, all SST studies suggest that the Pliocene/ late Miocene climate was characterized by a reduced zonal temperature gradient and a deeper/warmer thermocline relative to modern 1-4,6,9,10,31,32 -suggestive of a more El Niño-like background state or 'El Padre' 1, 5,12,33 . Any reduction in the zonal SST gradient would seem to contradict records of high primary production in the east Pacific at this time [13][14][15] .
While it is difficult to reconstruct nutrient delivery in upwelling areas, the boron isotopic (δ 11 B) composition of planktonic foraminifera provides a robust pH proxy for tracing changes in the chemical composition of seawater though time 34,35 . Critically, different paradigms of Pliocene climate, in which physical (SST) and biogeochemical (pH) gradients are either coupled or not, make distinct predictions about east and west tropical Pacific pH values and gradients (Fig. 1). Because surface pH values are strongly linked to regional and basin-scale ocean circulation patterns, which of these states the early Pliocene to late Miocene Pacific operated under has major implications for atmospheric circulation, weather patterns and ocean productivity.
The δ 11 B-pH proxy is based on the pH-dependent speciation of boron between two predominant aqueous species in seawater: boric acid (B(OH) 3 ) at low pH and borate ion (B(OH) 4 − ) at high pH, with an equilibrium fractionation of 10 B and 11 B between these two species. Only the borate ion is believed to be incorporated into biogenic carbonates and thus there is a predictable relationship between the δ 11 B of biogenic calcite (specifically, of planktonic foraminifera for the surface ocean) and ambient ocean pH 34,35 . Using boron isotopes, we here reconstruct zonal pH gradients in the tropical Pacific since the late Miocene (Fig. 2) and use a coupled ocean-atmosphere climate model to provide a mechanistic explanation for the apparent Pliocene paradox.
We analysed δ 11 B in at least ~150 tests per sample of the mixed layer-dwelling (see Extended Data Table 1) foraminifer Orbulina universa from the western equatorial Pacific (WEP, Ocean Drilling Program (ODP) Site 806) and eastern equatorial Pacific (EEP, ODP Site 846) at around 1 Ma, 3 Ma and 6 Ma. Each time slice includes four samples spanning about 300 kyr (n = 24). Our strategy was designed to sample a range of orbital configurations per time slice to capture the mean and variance in conditions. Although we do not expect to capture the full extent of orbital variability in each time slice, the variation in our ~3 Ma sample is comparable to that of an orbitally resolved record over the same time span, providing support for the approach (Fig. 2). We then converted δ 11  details of observational measurements of pH) to around 8.02 in the late Miocene (~6 Ma, mean pH uncertainty (2σ) ~0.05 from Monte Carlo simulations) (Fig. 2).
Our surface pH estimates for the WEP-a region understood to be approximately in equilibrium with the atmosphere 37 through the study interval due to a weak wind field and deep thermocline-are consistent with the decline one would expect in response to current estimates of atmospheric CO 2 concentrations during each interval (see Extended Data Fig. 1).
In the EEP (ODP Site 846), pH likewise declines from ~1 Ma to ~6 Ma but at a magnitude larger than the WEP, declining from an average modern pH value of around 8.007 to around 7.88 at about 6 Ma (uncertainty (2σ), ~0.06) (Fig. 2). This results in much greater zonal pH gradients (~2-3 times greater) at around 3-6 Ma than in recent times (average modern pH gradient ~0.06). The average zonal pH gradient (Fig. 2) varies from around 0.04 to around 0.17 to around 0.14 at approximately 1 Ma, 3 Ma and 6 Ma, respectively.
This three-fold increase in the equatorial Pacific zonal pH gradient at around 3 Ma and around 6 Ma is inconsistent with an El Niño-like coupling of physical (SST) and biogeochemical (pH) gradients (Fig. 1, coupled scenarios). Rather, our reconstructed zonal pH gradient suggests changes in zonal SST and pH gradients were decoupled over this interval (Fig. 1, decoupled). New simulations using a coupled oceanatmosphere climate model, meant to reproduce the background early Pliocene/late Miocene climate state 4,16 , run for the first time with active biogeochemical cycling (the Community Earth System Model, CESM 1.0.4, see Methods), generates pH changes similar to our results (Fig. 3, Extended Data Fig. 2a). In these simulations, cloud albedo has been modified to generate an early Pliocene/late Miocene climate state with strongly reduced zonal and meridional temperature gradients in agreement with SST reconstructions 4,10,16,[38][39][40][41] . While a simulation that reproduces early Pliocene/late Miocene climate without these cloud albedo modifications would be ideal, other (for example, Pli-oMIP2) models without these cloud albedo modifications target the mid-Pliocene warm period 26 . Future studies should therefore assess whether a simulation of early Pliocene/late Miocene climate without cloud albedo modifications and with more standard model parameters can replicate existing and our proxy data.
The modelled early Pliocene/late Miocene ocean has pH values of around 8.00-8.05 in the WEP and around 7.85-7.90 in the EEP (Fig. 3a) over a 25-55-m depth range (the approximate habitat of O. universa 36,42,43 ), similar to our δ 11 B-pH proxy data at ~6 Ma and, to a somewhat lesser extent, at ~3 Ma (Extended Data Fig. 2b-c). More broadly, however, model output over 25-75 m in depth generally agrees with δ 11 B-pH proxy data (Extended Data Fig. 2a, Extended Data Fig. 2d-g). A further comparison of the experiment-minus-control model anomalies and the proxy-minus-modern anomalies is made in Fig. 3b. Both data and model output show early Pliocene/late Miocene changes relative to the modern/ control of around −0.12 to -0.15 pH units in the east and relatively less pH change in the west. There is some model-data discrepancy in the magnitude of change in the west, with the model output exhibiting a slightly stronger anomaly of around −0.10 pH units relative to modern than is seen in the proxy data (around −0.05 pH units relative to modern; for further discussion see Supplementary Information). Nonetheless, the pattern in both the model output and proxy data is of increased zonal pH gradients at around 3 Ma and 6 Ma, even as zonal temperature gradients were decreased 1,3,7,8,11 (Figs. 2, 3; Extended Data Fig. 2).
The model simulations suggest that the enhanced ocean acidity in the east relative to the west can be explained by the upwelling of older waters in the east associated with a large-scale reorganization of Pacific circulation. The most dramatic changes include deep water formation in the North Pacific and an active Pacific meridional overturning circulation (PMOC) 16 , both of which are absent in the Pacific Ocean today and together increase the acidity of the intermediate waters being upwelled into the eastern equatorial Pacific. The modelled maximum PMOC strength is ~15 sverdrups (1 Sv = 10 6 m 3 s −1 , see Fig. 4b and Extended Data Fig. 3), comparable to the modern Atlantic meridional overturning circulation (AMOC) 44 .
Several features of this circulation regime contribute to the anomalously acidic water in the EEP cold tongue. First, deep water formation in the subarctic North Pacific introduces an anomalous interhemispheric gradient in mid-to deep-water pH relative to the pre-industrial ocean, with high-pH water in the Northern Hemisphere and low-pH water in the Southern Hemisphere (at the end of PMOC and AMOC flow paths) (Fig. 4a). Second, relatively more acidic (Fig. 4a) and older ( Fig. 4b) water is brought up by the PMOC to intermediate depths (1,000 m) in the southern and tropical Indo-Pacific (Fig. 4c). This acidic intermediate water is then readily upwelled in the EEP because of a warm, less stratified thermocline relative to modern 5,9 . According to Fig. 4, in the presence of the PMOC roughly 4 Sv of additional water enters the subtropical overturning cells (STCs) and equatorial thermocline, shoaling from around 1,000 m depth in the Southern Hemisphere to around 200 m in the tropics. This is close to 30% of the maximum volume transport of the PMOC. Third, along the pathways taken by water subducted in the subarctic North Pacific (Fig. 4d), water parcels spend a prolonged period of time at intermediate depths-taking ~800 years to arrive at the equator after transiting along the deep limb of the PMOC, with some parcels even taking several millennia-thereby gaining nutrients and acidity before upwelling in the EEP. In the modern climate, EEP thermocline waters largely originate from water subducting in the subtropical cells and only take around 10-60 years after subduction to transit into the EEP 45 .
While these results on their own cannot rule out potential effects of changes to other meridional overturning circulations (for example, the AMOC or shallow STCs) or water masses (for example, Antarctic Bottom Water, Antarctic Intermediate Water), a number of dynamical Article considerations make an active PMOC a robust interpretation of this pH record (see Supplementary Information). In brief, we expect that other overturning circulations and/or water masses are unlikely to drive the observed pH patterns owing to their depths (for example, ocean abyssal cells), locations of major impacts (the AMOC), or equilibration time scales and the sign of changes (for example, STCs) (Extended Data Fig. 3). The upwelling of old, acidic water in the east Pacific and associated delivery of more nutrients to the surface ocean contrasts with the view, based on modern observations, that nutrient delivery and productivity in warm oceans will be low, particularly in the eastern tropical Pacific, due to decreased wind-driven upwelling 46 . Our combined boron isotope data and model results (Fig. 3, Extended Data Fig. 2) show that old, more acidic, intermediate water can upwell in the EEP as a result of an active PMOC and decreased thermal stratification in the uppermost ocean. Our findings are most parsimoniously explained by more nutrient-rich water upwelling, and not greater upwelling volumes, since the STCs in the early Pliocene/late Miocene model run are clearly weaker than those in the pre-industrial control run (with northern and southern STCs of ~32 Sv and ~30 Sv respectively in the early Pliocene/late Miocene run versus ~38 Sv and ~46 Sv in the pre-industrial control run; see Extended Data Fig. 3e-h). Our results provide a dynamic mechanism for the enhanced biological productivity in the east Pacific despite warmer SSTs and a reduced zonal SST gradient during the early Pliocene to late Miocene, resolving the apparent paradox of decoupled SSTs and biological productivity.
The implications of our results for the future of the EEP depend on the timescale involved. Modelling results suggest such a more El Niño-like mean ocean state (in terms of SSTs) may take up to 100 years to develop in response to a rise in pCO 2 levels 21,47 . However, the long adjustment timescale of the deep ocean circulation implies that it would take on the order of ~1,000-1,500 years for the PMOC to develop in response to this warming (Extended Data Fig. 3c-d). In other words, the decoupling of temperature and pH gradients is predicted to develop progressively over a period of a thousand years or longer. With that in mind, increased temperatures and stratification in the EEP remain a major concern over the coming centuries given that EEP upwelling currently supports one of the most productive fisheries in the world and provides key economic and ecosystem services to local communities 48,49 , services whose importance will only increase with continued population growth and development. However, given that the emissions trajectories we set today will shape the climate for millennia 50 , such long-term dynamics as we've discussed here are not outside the realm of possibility either. More broadly, our results provide novel constraints on the structure and dynamics of the ocean of the early Pliocene and late Miocene by moving beyond the SST proxies typically employed and help constrain future changes in ocean dynamics as the ocean adjusts to a warming world.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-021-03884-7.   Our sampling strategy aimed to reconstruct the evolution of the zonal Pacific pH gradient since the latest Miocene. Samples were taken at three large time steps (~1 Ma, ~3 Ma and ~6 Ma) with four more closely spaced samples within each of these windows (every ~100 ka) to capture some extent of orbital variability. Our samples at ~1 Ma, for example, all fall within different stages of the glacial-interglacial cycles of the Pleistocene (Extended Data Fig. 4a, b). A more thorough treatment of glacial-interglacial pH changes is performed in other studies, particularly in Martínez-Botí and colleagues 64 . At a site in equilibrium with the atmosphere similar to our WEP site (PS2498-1 in the south central Atlantic), pH had a maximum range of 0.14 over the record (from ~2.24-15.90 ka) with ~95% of that variation (2σ) varying within just 0.06, for an average glacial-interglacial (~2.24-10.28 ka versus ~15.12-15.90 ka) pH change of 0.09. The EEP site ODP 1238 in ref. 64 has a maximum range and 2σ of 0.23 and 0.13 over the record (~1.03-20.60 ka). For comparison, our data at ~1 Ma (the time slice with the lowest within-site variance) has pH ranges of 0.05 and 0.06 for the WEP and EEP, respectively, suggesting-by comparison to ref. 64 -that we may not capture the full glacial-interglacial range in pH at 1 Ma. However, this does not affect the main finding of our study-that the zonal equatorial Pacific pH gradient in the early Pliocene/late Miocene (~0.17 and ~0.14 at ~3 Ma and ~6Ma, respectively) was enhanced relative to modern (~0.06)-given that the zonal gradients are roughly constant over the full glacial-interglacial cycle (for example, 0.08 in the glacial and 0.05 in the interglacial as calculated from ref. 64 ). Additionally, each of our samples are likely to integrate over ~10,000 years, minimizing concerns about mismatched time intervals captured by samples from the eastern and western Pacific sites (see also Extended Data Fig. 4).

Estimation of modern pH
To constrain modern values for eastern and western equatorial Pacific pH, we compiled observational measurements from all available cruise data near our sites from the past several decades to estimate the possible domain of modern variability in pH.
We have compiled direct measurements of pH (or occasionally, alkalinity and total dissolved inorganic carbon (DIC), from which pH may be calculated) from every cruise passing within 5° of latitude and longitude of our study sites, as available from the GLODAPv1 65 and v2 66 , PACIFICA 67 and WOCE 68 databases (18 cruises from 1992-2018, with 16 cruises in the west Pacific and two cruises in the east, not counting the eastern La Niña-year cruise, which was not taken near our study site (~20° west of it) so was only included for illustrative purposes and not included in any further analysis). Each cruise made multiple depth-profile measurements of either pH or alkalinity and DIC (68 profiles in total), which have been compiled into Extended Data Fig. 5. Where direct measurements of pH were not available, pH was calculated from alkalinity and DIC within Ocean Data View 69 using the Best Practices handbook 70 equilibrium constants, which give results in excellent agreement with other carbonate system calculations (for example, CO 2 SYS (ref. 71 )). The El Niño conditions under which each profile was taken was determined from the Oceanic Niño Index (ONI) product reported by NOAA's Climate Prediction Center 72,73 .
Estimates of modern pH in the eastern and western equatorial Pacific were derived only from neutral ENSO years (Extended Data Fig. 5a, with 12 profiles from one cruise in the east and 21 profiles from five cruises in the west, which were linearly interpolated onto uniform 1 m depth intervals and averaged together (Extended Data Fig. 5d). These average depth profiles were then averaged over the 25-50 m depth range (the approximate depth habitat of O. universa; see Supplementary Information and Extended Data Table 1) to produce estimates of modern pH against which to compare our proxy data (8.007 ± 0.062 and 8.067 ± 0.007 for the east and west, respectively (1σ)). An interesting observation from this pH compilation is that the zonal pH gradient is less pronounced at the surface than it is a few tens of metres deeper, making our choice to use O. universa more appropriate for detecting such signals in the past than, say, a shallower-dwelling species such as T. sacculifer.
More empirical data are needed from the eastern Pacific, and we therefore hesitate to put too much emphasis on data from a limited number of observations. For example, rather than the collapse of the zonal pH gradient expected during an El Niño event with its reduced upwelling, the pH gradient over 25-50 m depth that we find under El Niño conditions (~0.09, Extended Data Fig. 5b) is slightly greater but within variability of the average gradient found under neutral ENSO conditions (~0.06, Extended Data Fig. 5a). However, we do find an enhanced zonal pH gradient under La Niña conditions (~0.12, Extended Data Fig. 5c), consistent with increased east Pacific upwelling in a La Niña event (although note that the EEP cruise used to get this number was from a site ~20° west of our study site so is not directly comparable). Notably, we record in our proxy data zonal pH gradients even greater than this at ~3 Ma (~0.17 on average) and ~6 Ma (~0.14 on average).
Without more data available, the best practice is to constrain the possible domain of modern variability as we have done here. Accordingly, we have included in Fig. 2 both the average and 1σ of the 25-50 m averaged pH from all neutral-year profiles compiled from the east (blue diamond with error bars) and west (red diamond with error bars) Pacific. We have also reproduced Fig. 2 with every neutral-year profile's 25-50 m average pH being shown as dashes along its y axis (Extended Data Fig. 6).

Hypothesis figure
The schematic of Fig. 1 was generated by applying the same proportional reduction in the zonal thermal (SST) gradient (in the early Pliocene/late Miocene relative to modern) as found by different studies to the modern-day zonal pH gradient. First, modern WEP and EEP pH was found as described in the 'Estimation of modern pH' section, providing the anchor points on the left-hand side of the plot. Next, a representative early Pliocene/late Miocene WEP pH (~8.02) was estimated from an equilibrium state with 400 ppm atmospheric CO 2 and a modern-like alkalinity (2,275 μmol kg −1 ) (alkalinity is understood to have been relatively invariant over the Cenozoic, staying close to modern values 74 ). Finally, for the dotted and dashed blue lines, an early Pliocene/late Miocene EEP pH was estimated from the reported reductions in the late Miocene zonal thermal (SST) gradient relative to modern as found by Liu et al. 31 (dotted line, an ~82% reduction) and Zhang et al. 3 (dashed line, a ~28% reduction). The percent-reduction in the zonal thermal gradient relative to modern as found by each of those studies was then applied to the magnitude of the modern zonal pH gradient (~0.06), resulting in a representative early Pliocene/late Miocene zonal pH gradient of either 0.017 (dotted line) 31 or 0.049 (dashed line) 3 . Early Pliocene/late Miocene EEP pH was estimated by applying those pH gradients to the early Pliocene/late Miocene WEP pH estimate (~8.02), resulting in early Pliocene/late Miocene EEP pH values of ~8.003 (dotted) and ~7.971 (dashed). The solid blue line depicts the EEP average pH put out by the coupled ocean-atmosphere climate model used in this study in its early Pliocene/late Miocene-like run over the 25-55 m depth range (~7.886), which is interpreted to reflect a decoupling between early Pliocene/ late Miocene zonal thermal and pH gradients. The model could not be averaged down to 50 m exactly as its depth intervals went in steps of 10 m (that is, 25 m, 35 m, 45 m, 55 m).

Sample preparation for boron isotope and trace elements analysis
Initial preparation of marine sediment samples involved disaggregating samples in sodium hexametaphosphate solution (2 g l −1 ) buffered with ammonium hydroxide (NH 4 OH), washing in deionized water over a 64 μm sieve, and drying in a <50 °C oven. Samples were rewashed until visually clean. Individuals of O. universa were picked from the 300-355 μm size range. Occasionally individuals from other size fractions were included (down to 250 μm and up to 710 μm, but rarely outside of the 300-600 μm range) if foraminifera were sparse, as the boron-pH calibration for O. universa has no statistically significant size trend within this range 36 . Sample sizes typically range from ~150-200 foraminifera and ~2-4 mg of CaCO 3 .
All subsequent laboratory procedures were carried out in a Picotrace class ten clean lab at the Yale Metal Geochemistry Center (Yale University) in an over-pressurized hood equipped with boron-free filters and using Milli-Q water (18.2 MΩ with Q-Gard boron removal pack) and Teflon-distilled nitric acid (HNO 3 ). Approximately ~2.5 mg of solid carbonate sample (that is, the picked O. universa) was cleaned following the procedure outlined in Foster et al. 75 and Rae et al. 34 , and references therein. After lightly cracking open individual chambers between clean glass slides, the solid sample was repeatedly ultrasonicated and rinsed with Milli-Q water (5 × 30 s of ultrasonication) to remove clays. Organic matter was removed by subjecting the solid to oxidative cleaning using 250 μl of 1% hydrogen peroxide (H 2 O 2 ) in 0.1 M NH 4 OH at 80 °C for 3 × 5 min, with 15 s of ultrasonication in between. Afterwards, 250 μl of 0.0005 M HNO 3 was added to the solid for 30 s as a brief weak acid leach to remove any re-adsorbed contaminants. Samples were then rinsed twice with Milli-Q water. Finally, samples had 200 μl of Milli-Q water added to them before incrementally adding 0.5 M HNO 3 acid until the sample was totally dissolved (anywhere from 75 to 200 μl). The sample was centrifuged for 5 min to isolate any undissolved contaminant phases. A small aliquot (approximately 7% of the volume) of the supernatant was transferred to an acid-cleaned plastic centrifuge tube for trace metal analysis, and the remaining supernatant was transferred to a Teflon beaker for isotope analysis.
Boron was separated from the dissolved carbonate matrix by passing the sample through pre-cleaned heat-shrunk Teflon micro-columns containing 20 μl of ground, sieved (63-100 μm) boron-specific anionic exchange resin (Amberlite IRA 743 76 ) following the procedure of Foster and colleagues 75 . Before performing column chemistry, the columns were cleaned by passing a full-column's volume of 0.5 M HNO 3 plus 1 ml of 0.5 M HNO 3 plus 2 × 1 ml of Milli-Q through each of the columns. Samples were also buffered prior to column chemistry with boron-clean sodium acetate buffer (twice the volume of acid required to dissolve the solid sample) to maintain sample pH > 5. The sample matrix was then rinsed through the columns by the addition of 10 × 160 μl of Milli-Q water before being eluted with 5 × 120 μl of 0.5 M HNO 3 and collected in acid-cleaned Teflon beakers. Sample tails were also collected in separate Teflon beakers with a final elution of 1 × 120 μl of 0.5 M HNO 3 . Each batch of columns included a JCp-1 carbonate geostandard ( Japanese Geological Survey Porites Coral 77 ) to monitor column performance and reproducibility, as well as a total procedural blank (TPB) composed of 0.5 M HNO 3 and buffer to monitor cleanliness and potential sample contamination. The average blank contribution from the total procedural blank (TPB) included in each batch of columns was ~19 pg (~0.15% of average sample size).

Boron isotope analysis by MC-ICP-MS
Boron isotope analysis was performed at the Yale Metal Geochemistry Center on a Thermo Finnigan Neptune Plus MC-ICP-MS equipped with 10 12 Ω resistors and tuned before analysis to optimize maximum 11 B/ 10 B stability, following the procedure described of Foster et al. 75 and references therein. A Teflon barrel spray chamber and ammonia gas were used for boron washout 78 .
Prior to analysis, boron concentration and potential contamination by the Na-rich buffer were checked in tails an in 20 μl aliquots of samples diluted with 100 μl of 0.5 M HNO 3 . Tails typically made up <1% of the boron loaded. Depending on the sample size, boron isotope samples range from ~5-13 ng, yielding solutions for analysis of ~8-21 ppb (ng g −1 ).
Samples were bracketed with 50 ppb NIST SRM 951 boric acid standard, which was used to correct for machine-induced mass-fractionation and convert 11 B/ 10 B ratios to delta notation. To monitor accumulating blank contamination from the laboratory atmosphere while the vials were open for analysis, a blank consisting of the same volume of 0.5 M HNO 3 as the samples (600 μl) was analysed every third sample and used to correct sample δ 11 B values. Samples were measured in duplicate, with the mean δ 11 B being reported.
The average blank contribution from the total procedural blank (TPB) included in each batch of columns was ~19 pg (~0.15% of average sample size). Sample δ 11 B values were corrected using TPB δ 11 B values. The average blank correction (δ 11 B original - δ 11 B corrected ) was ~0.011‰ and therefore below typical measurement uncertainty.
The uncertainty of boron isotope measurements (2σ was typically 0.24‰) was determined from a relationship between signal size and external reproducibility: (1) 11 11 where [ 11 B] is the 11 B signal intensity in volts. This relationship is derived from repeat analysis of multiple batches of roughly mass-matched JCp-1 pass through columns with individual sample batches. For this study, the average δ 11 B value of JCp-1 was 24.18 ± 0.23‰ (2σ) (n = 8).

Trace element analysis by ICP-MS
Trace element analysis was performed at the Yale Metal Geochemistry Center on a Thermo Finnigan Element XR ICP-MS, again using a Teflon barrel spray chamber and ammonia gas and following the procedure outlined in Foster et al. 75 ).
pH and pCO 2 estimates from boron isotopes Calculation of pH from δ 11 B was performed using the O. universa calibration of Henehan et al. 36 , which shows good agreement with in situ pH in core tops from numerous geographically widespread sites covering a diverse range of oceanographic environments and ambient pH and δ 11 B of borate (Extended Data Fig. 4c). A Monte Carlo simulation was used to estimate uncertainty, following the example of numerous previous boron-pH studies 51,55,64,84,90 and using previously published code 35 .
Reported pH values are the mean of 10,000 replicates of a simulation in which sea surface temperature, sea surface salinity (SSS), δ 11 B calcite and δ 11 B SW , and alkalinity (for the purposes of calculating pCO 2 ) are randomly generated from within a frequency distribution spanning the uncertainty surrounding each parameter. Following the example of Chalk et al. 84 , we place normal distributions around SST, SSS, δ 11 B calcite and δ 11 B SW and a uniform distribution spanning a generous range (±200 μmol kg −1 ) around alkalinity. SSTs fed into the Monte Carlo simulations were derived from Mg/Ca data described in the 'Calculation of Mg/Ca SSTs' section and, unlike other proxy data, a wider and more conservative range of uncertainty of ±2 °C is applied instead of the root sum of squares of the analytical and calibration uncertainty (which was small, less than around 0.9 °C). Uncertainty around δ 11 B calcite is the reported analytical uncertainty (2σ). Because the residence time of boron in seawater is believed to be on the order of [10][11][12][13][14][15][16][17][18][19][20], the modern-day value of δ 11 B SW with its reported uncertainty (39.61 ± 0.04‰ (2σ)) is used here 94 . A generous range of sea surface salinity around a representative modern-day value is applied: 34.5 ± 1 psu. Reported uncertainty on pH is two standard deviations of the 10,000 simulated pH values. For the purpose of calculating pCO 2 estimates from the pH values, a modern-like value of alkalinity (2,275 μmol kg −1 ) was applied with a generous range of uncertainty (±200 μmol kg −1 ). Alkalinity is understood to have remained relatively constant over the Cenozoic 74 , making this a more robust assumption (but assuming a modern-like calcite saturation state gives similar results, Extended Data Fig. 1b). Finally, a map of modern pH (Extended Data Fig. 5e-g) is provided to show the locations of our study sites versus that of the other δ 11 B-derived pH data 51 in Fig. 2.

Physical and biogeochemical modelling
Numerical simulations were performed using the National Center for Atmospheric Research (NCAR)'s coupled ocean-atmosphere Community Earth System Model (CESM) version 1.2.2, with the T31 gx3v7 configuration designed for long-term paleoclimate simulations 95 . The atmospheric and land surface components (the Community Atmosphere Model 4 and Community Land Model 4) have a spectral truncation of T31, and the oceanic and sea ice components (POP2 and Community Ice Code) have a resolution ranging from 3° near the poles to 1° at the equator. Two simulations have been performed-a pre-industrial control simulation and an early Pliocene/late Miocene-like experiment that reproduces the magnitude and spatial structure of the large-scale warming patterns, specifically the reduction in the meridional and zonal gradients seen in Pliocene SST reconstructions 4,10,40,41 .
Both the control and early Pliocene/late Miocene-like experiment have the biogeochemical model enabled 96 , allowing us to evaluate the impact of the simulated changes in ocean temperature, salinity and circulation on ocean pH. The simulations have been run for 3,000 years, allowing the system to reach quasi-equilibrium (Extended Data Fig. 3a-d). The analysis presented here is based on the climatology of the last 100 years of each simulation (years 2901-3000), and Pliocene/ late Miocene anomalies are reported relative the pre-industrial control run.
The experimental design for the early Pliocene/late Miocene-like experiment is the same as the Early Pliocene-like simulation evaluated in Burls et al. 16 and Burls et al. 38,39 but with an active biogeochemical component. For full details, see Supplementary Information and Burls and colleagues 38,39 . In brief, mean cloud albedo is reduced in the extratropics and increased in the tropics by altering the liquid and ice water path in the shortwave radiation scheme in different latitudinal bands. The liquid water path of a cloud when it forms in the model is reduced poleward of 15°N and 15°S by 60%, while the ice and liquid water paths are increased equatorward of 15°N and 15°S by 240%. The cloud radiative forcing changes that arise due to the imposed water and ice path modifications result in a decrease in extratropical (8-90°N and 8-90°S) albedo of 0.04 and an increase in tropical albedo (8°S-8°N) of 0.06. This cloud albedo modification serves to compensate for under-represented cloud radiative effects from either poorly resolved cloud feedbacks (strongly positive in the extratropics and negative in the deep tropics; for example, refs. [97][98][99][100][101] or missing aerosol cloud indirect effects for example, Sagoo and Storelvmo 102 . All other boundary conditions (for example, the continental boundaries) and radiative forcings in the Pliocene/late Miocene-like experiment are the same as the pre-industrial control, that is, the default for the B_1850_BGC-BDRD CESM component set run on the T31_g37 grid 95 . The atmosphere, ocean, sea ice and land component of the climate system all freely evolve in response to the imposed changes in cloud shortwave radiative forcing reaching a near-equilibrium warming pattern that resembles early Pliocene and late Miocene SSTs. Note that atmospheric CO 2 was left at the pre-industrial value of 284.7 ppm in the atmospheric component and the global warming seen in this simulation is instead supported by the cloud albedo changes and the associated feedbacks, for example, the water vapour feedback (in reality elevated Pliocene CO 2 concentrations were likely to be the primary forcing responsible for the Pliocene warmth with the cloud radiative forcing changes resulting as a feedback, or due to different Pliocene aerosol concentrations). Atmospheric CO 2 was set to 400 ppm in the ocean biogeochemistry (BGC) component to correctly simulate the implications of elevated Pliocene CO 2 concentrations on ocean biogeochemistry and pH.

Lagrangian analysis of water parcel trajectories
In this study we used the Lagrangian analysis of water parcel trajectories of Thomas and colleagues 103 . Water parcel trajectories were modelled to identify the subsurface pathways taken by the PMOC water. This involved identifying both the surface source regions of deep-water formation in the North Pacific and the mixed layer destinations where that water resurfaces, and tracing the pathways between them, determined according to a Lagrangian ocean analysis approach that tracks water parcels in the Pacific Ocean of the CESM climate model (see the 'Physical and biogeochemical modelling' sectionabove for model description). This analysis was performed using the Ariane software package 104 , which models water parcel trajectories by time-integrating velocity fields output by the ocean model. Trajectories calculated from ocean model momentum equations in this way-when coupled with a parameterization for subgrid-scale eddy-induced velocities (EIVs) 105 -depict the purely advective pathways of water parcels/tracers (that is, without the effect of diffusion). To resolve the long residence timescales of PMOC water in the deep ocean, particles were allowed to loop continuously over the 100 years of model output (see van Sebille et al. 106 for a discussion on particle looping).
Model tracer particles were seeded in every grid cell of the Pacific Ocean basin at 25°N and at every timestep (~1 million particles) before being integrated both backwards in time (to identify the particle mixed layer source regions at high northern latitudes) and forwards in time (to identify their mixed layer destinations). Time-integration in both cases was stopped when particles reached the model mixed layer. Particles found to originate close to the northern margin of the Pacific Ocean basin (north of ~40°N), clearly identifiable as PMOC water that is separate from the subtropical sources of the shallow wind-driven cells, were then isolated and subsequently run forwards in time from their original starting positions at 25°N. Of these, the particles that upwelled into the tropical Pacific (approximately 15% of them) were then further identified and their pathways and timescales analysed. A random subset of 12,000 particles (a number manageable within computer memory limitations) was then selected to identify the pathways taken between the sources and tropical sinks of PMOC water (Fig. 4c).
Because the CESM ocean model used in this study is configured on an Arakawa B-grid, ocean model velocities had to first be linearly interpolated onto a C-grid configuration 107 to be compatible with the Ariane software. EIV values also had to be set to zero at the ocean boundary to prevent unrealistic propagation of particles; however, the Pacific Ocean domain was found to be minimally affected by this problem, and results overall were only weakly sensitive to it.

Data availability
The proxy data and model output produced in this study are available as .xlsm and .nc files in NOAA's paleoclimate data repository (https:// www.ncdc.noaa.gov/paleo/study/33252) (https://doi.org/10.25921/ AMPV-J413). Source data are provided with this paper.