Decoupling of temperature and pH gradients along the equatorial Pacific during the Pliocene

Decoupling of temperature and pH gradients along the equatorial Pacific during the Pliocene Authors: Madison G. Shankle , Natalie J. Burls, Alexey V. Fedorov, Matthew D. Thomas, Donald E. Penman , Heather L. Ford, Peter Jacobs, Noah J. Planavsky, Pincelli M. Hull Affiliations: Department of Earth and Planetary Sciences, Yale University, New Haven, CT, 06511, USA. 5 Department of Atmospheric, Oceanic & Earth Sciences, George Mason University, Fairfax, VA, 22030, USA. University Corporation for Atmospheric Research/GFDL, Princeton, NJ, 08540, USA. School of Geography, Queen Mary University of London, London, E1 4NS, UK. Department of Environmental Science and Policy, George Mason University, Fairfax, VA, 22030, 10 USA. *Correspondence to: mgs23@st-andrews.ac.uk, pincelli.hull@yale.edu. †Current address: School of Earth & Environmental Sciences, University of St Andrews, St Andrews, KY16 9AL, UK. ††Current address: Department of Geosciences, Utah State University, Logan, UT, 94321 USA. 15

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 Pliocene Pacific operated under has major implications for atmospheric 95 circulation, weather patterns, and ocean productivity. The boron-pH proxy therefore offers a potent tool for constraining both Pliocene tropical climate and ocean dynamics.
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 B 10 and B 11 between these two species. Only 100 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 29,30 . Using boron isotopes, we thus reconstruct zonal pH gradients in the tropical Pliocene Pacific (Fig. 2) and use global circulation models to provide a mechanistic explanation for the apparent Pliocene paradox. 105 Here we analyzed  11 B in ~100-shells per sample of the near surface-dwelling foraminifer Our pH values give atmospheric pCO2 estimates of 331 ± 48 ppm [2] at 1Ma, 434 ± 64 ppm at 3Ma, 439 ± 64 ppm at 6Ma when assuming a generous range around the approximate modern-day value for alkalinity (2275 ± 200 mol/kg 32 ) (Fig. S1). Alkalinity is understood to have remained relatively invariant over the Cenozoic 33 , making it a robust assumption to use, however the assumption of a modern-like value for the calcite saturation state with a generous range ( = 5 ± 120 2) gives similar results (Fig. S2).
This three-fold increase in the Pliocene equatorial Pacific zonal pH gradient 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 130 and pH gradients were decoupled in the Pliocene (Fig. 1, decoupled). New simulations using a coupled ocean-atmosphere climate model with active biogeochemical cycles (the Community Earth System Model, CESM 1.0.4, see Methods), which reproduces large-scale surface temperature gradients resembling Pliocene SST reconstructions 16,34,35 , generate pH changes similar to our results (Fig. 3). The modeled Pliocene ocean has pH values of ~8.00-8.05 in the 135 WEP and ~7.80-7.85 in the EEP (Fig. 3a) at 55m depth (the approximate habitat of O. universa [36][37][38] ), similar to our δ 11 B-pH proxy data. Model output over 25-75m depth generally agrees with δ 11 B-pH proxy data, however (Figs. S3, S4). This proxy-model agreement provides compelling evidence for increased zonal pH gradients in the Pliocene even as zonal temperature gradients were decreased 1,5,6,9,11 (Figs. 2, 3; Figs. S4, S5). 140 The Pliocene model simulations suggest that the enhanced ocean acidity in the east relative to the west can be explained by a large-scale reorganization of Pacific circulation associated with 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. The modeled Pliocene PMOC strength is ~17 Sverdrups (1 Sv = 10 6 m 3 s -1 ), comparable to the modern Atlantic meridional 145 overturning circulation (AMOC). Several features of the Pliocene PMOC contribute to the anomalously acidic water in the EEP cold tongue. First, deep water formation in the North Pacific introduces an anomalous interhemispheric gradient in mid-to deep water pH relative to the preindustrial 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, PMOC 150 brings deeper, low-pH water to intermediate depths (1000 m) (Fig. 4b). This acidic intermediate water is then readily upwelled in the EEP because of a warm, less stratified thermocline relative to modern 3,7 . A Lagrangian analysis of the trajectories of water parcels subducted during North decreased wind-driven upwellin 17 . Our coupled boron isotope data and model results show that old, more acidic, intermediate water can upwell in the EEP as a result of PMOC and decreased thermal stratification in the uppermost ocean. Our results provide a dynamic mechanism for the enhanced biological productivity in the east Pacific despite warmer SSTs and a reduced zonal SST 165 gradient during the Pliocene, resolving the apparent paradox of decoupled SSTs and biological productivity.
The implications of our results for the future of the EEP depend on the time scale involved.
Modeling results suggest such an El Niño-like mean ocean state may take roughly ~100-200 years to develop in response to a rise in pCO2 levels 27,40 . However, the long adjustment timescale of the 170 deep ocean circulation implies that it would take on the order of ~1000-1500 years for PMOC to develop in response to this warming (Fig. S4). In other words, the decoupling of temperature and pH gradients is predicted to develop progressively over a period of a thousand years or longer.
This warming and stratification in the EEP is still of concern since EEP upwelling currently supports one of the most productive fisheries in the world and provides key economic and 175 ecosystem services to local communities 19,20 . More broadly, our results provide novel constraints on the structure and dynamics of the Pliocene ocean by moving beyond the SST proxies typically employed and help constrain future changes in ocean dynamics and the oceans adjust to a warming world.

Study Sites and Species Selection
To reconstruct the zonal pH gradient across the equatorial Pacific, samples were taken at where waters are well mixed and in communication with the atmosphere, and has been calibrated for boron isotopes 37 , making it a suitable candidate for this study of mixed layer pH 36-38 . Figure  495 The schematic of Fig. 1 was generated by applying the same proportional reduction in the zonal thermal (SST) gradient (in the Pliocene relative to modern) as found by different studies to the modern-day zonal pH gradient. First, pre-Industrial pH for both the WEP and EEP was taken in the zonal thermal gradient relative to modern as found by each of those studies was then applied to the magnitude of the pre-Industrial zonal pH gradient (~0.5) 31 , resulting in a Pliocene zonal pH gradient of either 0.014 (dotted line) 41 or 0.041 (dashed line) 11 . Pliocene EEP pH was estimated by applying those Pliocene pH gradients to the Pliocene WEP pH estimate (~8.02), resulting in Pliocene EEP pH values of ~8.006 (dotted) and ~7.979 (dashed). The solid blue line depicts the 510 EEP pH put out by the earth system model used in this study in its Pliocene-like run (~7.82), which is interpreted to reflect a decoupling between Pliocene zonal thermal and pH gradients.

Calculation of Mg/Ca SSTs
Planktonic foraminiferal Mg/Ca was measured on an aliquot of each δ 11 B sample and used to estimate SSTs in order to accurately calculate pH from δ 11 B. Given the residence times of magnesium (~13 Myr) 57  Carlo simulation in which the involved parameters (sea surface temperature, sea surface salinity (SSS),  11 Bcalcite, and  11 BSW, and alkalinity (for the purposes of calculating pCO2)) are randomly generated from within either a normal or uniform frequency distribution spanning the uncertainty surrounding each parameter. Following the example of Chalk et al., 2017, 62 , we place normal 620 distributions around SST, SSS,  11 Bcalcite, and  11 BSW and a uniform distribution spanning a generous range (±200 mmol/kg) 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 o C is applied instead of the root sum of squares of the analytical and 625 calibration uncertainty (which was small, <~0.9 o C). Uncertainty around  11 Bcalcite is the reported analytical uncertainty (2). Because the residence time of boron in seawater is believed to be on the order of 10-20 Myr 70-72 , the modern-day value of δ 11 BSW with its reported uncertainty (39.61 ± 0.04‰ (2)) is used here 42 . A generous range of sea surface salinity around a representative modern-day value is applied: 34.5 ± 1 psu. Reported uncertainty on pH is one standard deviation 630 of the 10,000 simulated pH values. For the purpose of calculating pCO2 estimates from the pH values, a modern-like value of alkalinity (2275 mmol/kg) was applied with a generous range of uncertainty (±200 mmol/kg). Alkalinity is understood to have remained relatively constant over the Cenozoic 33 , making this a more robust assumption than, for instance, the assumption of a modern-like calcite saturation state (although this gives similar results, Fig. S2). 635

Effect of SST Estimates on pH Reconstructions
To calculate pH, SSTs were estimated from planktonic foraminiferal calcite Mg/Ca using pH values produced by different SST proxies increases going back in time at SST records diverge more and more but still rarely surpasses ~0.05 in our data (Fig. S8). Similarly, calculating pH using the SSTs derived from the BAYMAG calibration 67 (Fig. S10)-while shifting pH values down 660 ~0.05 at the 3Ma and 6Ma time points-does not affect the main conclusions of this study. conclusions of this paper (compare Fig. 4 and Fig. S11).

Lagrangian Analysis of Water Parcel Trajectories
Water parcel trajectories were modeled in order to identify both the surface source regions of deep water formation in the North Pacific and the mixed layer destinations where that water resurfaces, determined according to a Lagrangian ocean analysis approach that tracks water parcels 710 in the Pacific Ocean of the CESM climate model (see above for model description). This analysis was performed using the Ariane software package 78  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 81 in order 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 735
The data generated and analyzed in this study are included in this published article and in its supplementary Extended Data section.

Acknowledgements:
The authors acknowledge and appreciate the aid of Jamie Robbins, Wayne Strojie, Dan 740 Asael, and Shuang Zhang in supervising clean lab chemistry, boron and trace element analysis, and data processing. The authors also appreciate discussion with James Rae concerning the Monte  pCO2 estimates (blue markers and error bars) made using pH proxy data of this study from the 805 WEP, a region understood to be in equilibrium with the atmosphere 25 , assuming a modern-like alkalinity of 2275 mmol/kg. Error bars depict uncertainty (2) returned from a Monte Carlo simulation (n = 10,000 runs) including uncertainty from the pH estimate (i.e. uncertainty on measured  11 B, the  11 B composition of seawater 42 , SST, and sea surface salinity (SSS)), and ± 200 mmol/kg on alkalinity. Grey data points show  11 B-deried estimates of pCO2 from previously 810 published studies, with shading showing reported confidence intervals 43,62,82-86 .   pH as derived from the model output (lines and shaded areas) and from  11 B measurements (circle markers, error bars show uncertainty (2) returned by a Monte Carlo simulation, as in Figure 2B).
Markers are placed at their approximate longitudes (latitudes are 0 o N for western site ODP 806 and -3 o N for eastern site ODP 846). O. universa is understood to be a mixed-layer species 36-38 with 840 a depth habitat ranging from ~0-75m, although it has been found at deeper depths 87,88 . The pre-Industrial longitudinal gradient in surface pH is included for reference 31 .