Coupling of Surface Ocean Heat and Carbon Perturbations over the Subtropical Cells under Twenty-First Century Climate Change

Itiswellestablishedthattheoceanplaysanimportantroleinabsorbinganthropogeniccarbon C ant fromthe atmosphere. Under global warming, Earth system model simulations and theoretical arguments indicate that the capacity of the ocean to absorb C ant will be reduced, with this constituting a positive carbon–climate feedback. Here we apply a suite of sensitivity simulations with a comprehensive Earth system model to demonstratethatthesurfacewatersoftheshallowoverturningstructures(spanning45 8 S–45 8 N)sustainnearly half of the global ocean carbon–climate feedback. The main results reveal a feedback that is initially triggered by warming but that ampliﬁes over time as C ant invasion enhances the sensitivity of surface p CO 2 to further warming, particularly in the warmer season. Importantly, this ‘‘heat–carbon feedback’’ mechanism is distinct from (and signiﬁcantly weaker than) what one would expect from temperature-controlled solubility perturbations to p CO 2 alone. It ﬁnds independent conﬁrmation in an additional perturbation experiment with the same Earth system model. There mechanism denial is applied by disallowing the secular trend in the physical state of the ocean under climate change, while simultaneously allowing the effects of heating to impact sea surface p CO 2 and thereby CO 2 uptake. Reemergence of C ant along the equator within the shallow over-turning circulation plays an important role in the heat–carbon feedback, with the decadal renewal time scale for thermocline waters modulating the feedback response. The results here for 45 8 S–45 8 N stand in contrast to what is found in the high latitudes, where a clear signature of a broader range of driving mechanisms is present.


Introduction
An important priority in climate research is to identify and understand processes that modulate or limit the rate at which the ocean can absorb anthropogenic carbon from the atmosphere under future climate change.Recent analyses demonstrate that over 1959-2017, approximately 25% of anthropogenic carbon emissions were absorbed by the ocean (Friedlingstein et al. 2019), but perturbations to the ocean state in the future may limit the ocean's ability to absorb anthropogenic carbon.For the case where anthropogenic perturbations to the ocean state may result in more CO 2 being left in the atmosphere and a resulting radiative heating of the Earth system, we refer to this as a positive climate feedback.There is long-standing interest in how anthropogenic heat entering the ocean can trigger climate feedbacks through reduced ocean CO 2 uptake (Banks and Gregory 2006;Xie and Vallis 2012;Frölicher et al. 2015;Winton et al. 2013).To date, these studies have emphasized potential for climate feedbacks through the impact of perturbations on the physical state of the ocean that impact ventilation of heat and carbon anomalies differently.
With the availability of Earth system models, including prognostic ocean biogeochemistry components, a series of related efforts over the past 15 years have focused on developing a framework for quantifying carbon-climate feedbacks.Efforts to identify carbon-climate feedbacks with climate models began with the studies of Maier-Reimer et al. (1996), Sarmiento et al. (1998), and Joos et al. (1999), where the focus was on marine carbon-climate feedbacks.This was subsequently considered by Cramer et al. (2001) for land carbon cycle-climate feedbacks.A now ubiquitous framework for quantifying feedbacks was developed in a series of studies that addressed both marine and land contributions, first by Cox et al. (2000) and then by Friedlingstein et al. (2003Friedlingstein et al. ( , 2006)).Although the original studies focused on global diagnostics of feedbacks, subsequent work (an example being the multimodel comparison study of Roy et al 2011) considered a regional breakdown of where ocean carbon-climate feedbacks find strongest expression.By considering four models that were forced under a Special Report on Emissions Scenarios (SRES) A2 future  scenario (Nakic ´enovic ´and Swart 2000), the study of Roy et al. (2011) indicated that the low-latitude regions spanned by 448S-498N dominate global carbonclimate feedbacks by contributing 71% of the global total when averaged across the models.This is important as it includes neither of the high-latitude regions that were emphasized in the pioneering works of Maier-Reimer et al. (1996) (North Atlantic focus) and Sarmiento et al. (1998) (Southern Ocean focus).
Ocean carbon feedbacks on the climate system under anthropogenic carbon emissions have been addressed historically in two different ways, both addressing how ocean perturbations may result in more CO 2 in the atmosphere than would be expected for an unperturbed ocean state.The first and early framing relates to thermodynamic perturbations to carbon chemistry and the potential impact of CO 2 invasion on the CO 2 buffering capacity of surface waters (Revelle and Suess 1957;Zeebe and Wolf-Gladrow 2001).As the CO 2 buffering capacity decreases with increasing CO 2 invasion, the ocean's ability to absorb CO 2 from the atmosphere relative to an unperturbed ocean state will be diminished, and such feedbacks could in principal occur for an unperturbed (preindustrial) circulation state of the ocean.
The second class of feedback mechanisms, developed over subsequent decades, invoke perturbations to the physical circulation state of the ocean, and ensuing perturbations to mixing, biology, and a range of other consequences.If the upper ocean warms and the surface polar regions freshen, increased stratification reduces the ocean interior ventilation rate.Together with a reduced solubility of CO 2 in warmer waters, it is expected that this would reduce the uptake of CO 2 (Siegenthaler and Oeschger 1978;Maier-Reimer et al. 1996;Sarmiento et al. 1998;Joos et al. 1999;Schwinger and Tjiputra 2018).Within the marine carbon cycle feedback framework that has gained widespread usage in the climate literature (Cox et al. 2000;Friedlingstein et al. 2003Friedlingstein et al. , 2006;;Roy et al. 2011), Earth system models have been applied to quantify both globally and regionally the degree to which ocean uptake of CO 2 is decreased for a perturbed relative to an unperturbed preindustrial state of the ocean.Although the diagnostics have proven invaluable for the quantification of feedbacks, there is not yet a consensus on the relative contributions of physical drivers (subduction and ventilation), perturbations to biology, and drivers that operate through perturbations to the carbon dioxide buffering capacity of surface waters.
The main objective of this study is to demonstrate that a heat-carbon feedback over the low latitudes (458S-458N) is triggered through heating of surface waters, with amplification through time occurring through the cumulative invasion flux of anthropogenic carbon boosting the sensitivity of sea surface pCO 2 to warming, with this amplification being larger in the warmer season.This is investigated through a suite of perturbation sensitivity studies with a CMIP5-generation Earth system model, using a high greenhouse gas emissions scenario (historical/RCP8.5)as the baseline for evaluating changes over the period 1861-2100.Although our principal interest is in mechanisms occurring over the low latitudes, we present global diagnostics in order to explore potentially contrasting mechanisms that occur over low and high latitudes in sustaining marine carbon-climate feedbacks.

a. Simulations
Here we apply a suite of five simulations with the Geophysical Fluid Dynamics Laboratory (GFDL) Earth system model (ESM2M; Dunne et al. 2012Dunne et al. , 2013) ) under historical and future atmospheric CO 2 concentration pathways.The physical state model underlying ESM2M is an updated version of GFDL's coupled model CM2.1 (Delworth et al. 2006).The ocean component is version 4.1 of the Modular Ocean Model (MOM4p1) (Griffies 2009) with approximately 18 horizontal resolution in the ocean that is enhanced meridionally near the equator and with 50 vertical layers.The atmospheric component of ESM2M is AM2 (Anderson et al. 2004), with approximately 28 horizontal resolution in the atmosphere.The ocean biogeochemical model is Tracers of Ocean Phytoplankton and Allometric Zooplankton code, version 2 (TOPAZv2), including 30 tracers to represent cycles of carbon, oxygen, the major macronutrients, and iron (Dunne et al. 2010).The fidelity of ESM2M to observational constraints over the historical period has been addressed by Dunne et al. (2012Dunne et al. ( , 2013)), and the small transient climate response (TCR) of ESM2M relative to the CMIP5 models has been described by Forster et al. (2013).The fidelity of ESM2M in simulating ocean uptake of anthropogenic carbon has been described in the study of Frölicher et al. (2015).
The experimental design underlying this study is intended to deconvolve the relative contributions of thermodynamics of marine carbon chemistry and ventilation rate changes (ocean physics) and biological perturbations in sustaining marine carbon-climate feedbacks.Five simulations are organized to interpret feedbacks on the fully coupled simulation (COU) that follows a historical/RCP8.5 pathway over 1861-2100(van Vuuren et al. 2011) (Fig. 1).The simulations apply prescribed mixing ratios for atmospheric CO 2 and other non-CO 2 radiative forcing agents.The first principal simulation (COU) uses the historical/RCP8.5CO 2 boundary conditions for radiative as well as for the land and ocean biogeochemistry components.The second simulation [biogeochemically coupled (BGC)] is the same as the COU simulation, except the concentrations of greenhouse gases and other radiative agents were kept at preindustrial levels in the radiation code, while atmospheric CO 2 concentrations follow the historical/ RCP8.5 pathway in the gas exchange code.Conversely, the third simulation [radiatively coupled (RAD)] includes the increase of greenhouse gases and their radiative agents on the physical state of the climate system with associated perturbations to the global carbon cycle, but for air-sea gas exchange the ocean sees a persistent preindustrial mixing ratio for atmospheric CO 2 .We also included a control simulation under preindustrial conditions (PIN), with this consisting of a 240-yr extension of the original spinup run.
An additional simulation, COU constant circulation (COU_CC), was performed using precisely the same physical state evolution as was used for BGC over the corresponding period 1861-2100 period (a preindustrial ocean state), but the thermodynamic terms impacting sea surface pCO 2 used for gas exchange were calculated using high-frequency output from the COU run for which there are important climate transients in surface temperature and salinity.This was accomplished by first saving daily mean SST and SSS from the COU simulation over the full interval 1861-2100.These daily mean SST and SSS fields were then read into the gas exchange routines for the calculation of pCO 2 in the COU_CC simulation, and linearly interpolated to the (2 h) time step of the ocean model, thereby impacting gas exchange at the sea surface through their impact on DpCO 2 .The other terms required in the calculation of pCO 2 at the surface in the COU_CC simulation, namely surface dissolved inorganic carbon, surface alkalinity, PO 4 , and SiO 4 , were taken directly from the fully prognostic fields in the biogeochemistry model (TOPAZ) in the COU_CC simulation.
In this way, the CO 2 system in surface seawater (including the buffering capacity) for the new COU_CC run feels the warming transient, but COU_CC is not impacted by the changes in stratification and ventilation rates/pathways experienced by the COU run itself.Specifically in the CO 2 chemistry calculations (Mehrbach et al. 1973;Najjar and Orr 1999), K 0 (the solubility constant) is perturbed through the effects of temperature and FIG. 1. Schematic for five model runs considered here, following the nomenclature of Gregory et al. (2009).The fully coupled case is denoted COU, the radiatively coupled case is denoted RAD, the biogeochemically coupled case is denoted BGC, and the preindustrial run is denoted PIN.An addition run is introduced in our analysis here, namely COU_CC (COU with constant climate, namely a preindustrial circulation state).For COU_CC the circulation state of the ocean is that shared by the PIN and BGC, but data override with the daily SST and SSS fields from the COU run used to calculate pCO 2 for COU_CC within the gas exchange routines.As such, COU_CC should be interpreted as being analogous to the COU case except with mechanism denial in terms of not allowing a transient in the circulation state.salinity, as are K 1 and K 2 (dissociation constants of carbonic acid).Thus, the COU_CC perturbation experiment is best understood as being a COU analog, where through mechanism denial the transient in the ventilation rate for the ocean has been disallowed.It is important to keep in mind that the cumulative effect of the perturbation applied in COU_CC results in reduced global carbon uptake relative to BGC (the unperturbed case).It is also worth noting that we chose to build COU_CC on the preindustrial circulation state of BGC with targeted perturbation to SST and SSS in the gas exchange routines, rather than the reverse configuration building on a transient circulation state with targeted SST and SSS perturbations from a preindustrial state of the model.Our choice was motivated by the view that the results based on a preindustrial circulation state would be much more straightforward to interpret.

b. Data-based sea surface pCO 2 products
Our primary analysis here focuses on model projections, but we also make use of observational products to assess the fidelity of ESM2M in several regions of interest.We utilize two databased pCO 2 products: The Princeton University Markov chain Monte Carlo (PU-MCMC) method (Majkut et al. 2014) and the product of the Japan Meteorological Agency (JMA) (Iida et al. 2015).To represent the evolution of pCO 2 in SSTsalinity-normalized dissolved inorganic carbon (nDIC) space for the observational products PU-MCMC and JMA, it is necessary to specify a second variable in the carbonate system of seawater to derive DIC.We have chosen here to use two published empirical methods for calculating alkalinity to this end.For the two Pacific sites [the formation region for North Pacific Subtropical Mode Water (NPSTMW) and the equatorial Pacific divergence region (EQDIV)], we use the interannually varying monthly climatology of Takatani et al. (2014), and for the subpolar North Atlantic we use the product of Lee et al. (2006) [the method of Takatani et al. (2014) has not yet been extended to include the Atlantic].As was presented in the study of Takatani et al. (2014), incorporating remotely sensed sea surface height variations in an empirically based account of surface alkalinity offers a number of distinct advantages over the method of Lee et al. (2006), in particular in regions such as NPSTMW that are on the poleward flank of the subtropical gyres.To maintain consistency with the observational products used to derive the product of Takatani et al. (2014), we also use sea surface salinity from the product of Usui et al. (2006) and the sea surface temperature product of Kurihara et al. (2006) in calculating nDIC over the 2000s, and subsequently for the construction of a climatology over the same period.
Likewise, the same temperature and salinity fields are provided for the North Atlantic to the calculation of alkalinity using the method of Lee et al. (2006).To facilitate comparison, the figures normalize all fields to constant alkalinity and salinity values that are given in the figure panels (this will be shown in Figs. 5 and 6).All of the CO 2 chemistry calculations are performed using the method of Lueker et al. (2000), as are the isolines of constant pCO 2 in each of the panels.c.Diagnostics of the changes in the CO 2 buffering capacity of surface waters Additionally, we consider the Revelle factor (Revelle and Suess 1957) as a diagnostic of the buffering capacity of CO 2 in surface seawater, with this commonly expressed with Revelle factor 5 (DpCO 2 /pCO 2 )/(DDIC/DIC) .
For the calculations used here, we employ the CO2SYS algorithms of van Heuven et al. (2011), to which we first calculate the monthly mean Revelle factor using DIC concentrations, in conjunction with SST and SSS (pCO 2 is also required for the iterative calculations performed offline with CO2SYS).

Results
We begin by considering the integrated air-sea CO 2 fluxes over the period 1861-2100 for the different simulations and regions (Fig. 2).The global ocean CO 2 uptake for the COU case (black solid line in Fig. 2a) is reduced by 1.3 PgC yr 21 during the 2090s in comparison with the BGC case (green line in Fig. 2a), indicating the well-known positive global carbon-climate feedback, with the contrast between COU and BGC being larger than the contrast between RAD (red line in Fig. 2a) and PIN (blue line in Fig. 2a) of 0.19 PgC yr 21 .This nonlinearity has previously been discussed in the studies of Gregory et al. (2009) and Schwinger et al. (2014).
We next consider integrated CO 2 uptake over three distinct latitude bands of the global domain.Given our primary interest in CO 2 uptake over the shallow overturning cells (encompassing thermocline and subpolar mode waters) spanning 458S-458N (Iudicone et al. 2016;Toyama et al. 2017), we divide the global domain into three bands.These are 458S-458N (shallow overturning), 458-908N (northern high latitudes), and 908-458S (southern high latitudes).A comparison over the three latitude bands reveals that 60% of the global cumulative carbon uptake occurs within the low latitudes (458S-458N) by the end of the twenty-first century (Figs.2b,c).In fact, this dominance of thermocline density-class waters in uptake is not surprising considering that the latitudes 458S-458N span 82% of the global ocean surface area.The regions spanning 908-458S and 458-908N absorb approximately 30% and 10% of the global total, respectively.
Importantly, the spatially integrated CO 2 fluxes for COU_CC over 458S-458N (dashed line in Fig. 2c) reveals that surface warming in conjunction with the invasion flux of CO 2 is able to drive nearly all of the carbon-climate feedback over this region.Thus, the large-scale reduction in CO 2 uptake over 458S-458N for COU relative to BGC is first triggered by surface heating, and then amplified in time by the invasion flux of anthropogenic CO 2 boosting the sensitivity of surface temperature perturbations.In other words, the underlying driver of changes in pCO 2 are SST rather than DIC or alkalinity (ALK) driven.This stands in contrast to the higher latitudes of both hemispheres (Fig. 2b), where perturbations to the ocean circulation state under climate change sustain important changes to the natural carbon cycle and clearly to pCO 2 (DIC-driven changes) (de Lavergne et al. 2014;Winton et al. 2013).
The reduction in cumulative anthropogenic carbon uptake due to climate change (differences between the COU and BGC simulations) over 1861-2100 is investigated next.Cumulative carbon uptake is reduced by 10% (567 PgC in COU vs 632 PgC in BGC, with a total perturbation of 64.2 PgC) due to climate change by 2100 (black line in Fig. 2d).The low latitudes contribute 48% (31 PgC) to the global difference, whereas 908-458S and 458-908N contribute only 38% (24 PgC) and 14% (9 PgC), respectively (Fig. 2d).Although the contribution of the Southern Ocean is larger per unit area than the latitudes spanned by the shallow overturning, the vast expanse of the region spanning 458S-458N allows this region to contribute nearly half of the global total.Over global scales, it warrants mention that by 2095, the large ensemble suite of 30 runs with ESM2M (Rodgers et al. 2015), identical in their forcing and model configuration to the COU run here, differ in their global total DIC inventories with a standard deviation of 1.2 PgC, with this being only 2% of the 64.2 PgC difference in global total DIC inventories for BGC and COU.This confirms that the difference in total DIC inventories between BGC and COU unambiguously represents an emerged response to the applied perturbation, rather than reflecting differences associated with natural variability.FIG. 2. For ESM2M over 1861-2100, (a) globally integrated uptake for COU, RAD, BGC, PIN, and COU_CC cases as defined in Fig. 1b combined uptake integrated over both (b) the high northern latitudes 458S-458N and the Southern Ocean defined as the region 908-458S, and (c) uptake integrated over 458S-458N.(d) The differences in the cumulative uptake for BGC and COU globally and by latitude range.
To gain insight into the mechanistic controls on carbon-climate feedbacks for our three regions of interest, we turn our attention to the annual mean perturbations for COU relative to BGC (considered as averages over 2080-99) for SST (Fig. 3a), nDIC (Fig. 3b), pCO 2 (Fig. 3c), CO 2 fluxes (Fig. 3d), changes in the surface ocean Revelle factor (Fig. 3e), and finally changes in surface ocean DIC concentrations (Fig. 3f).SST perturbations in ESM2M are characterized by broad warming over most of 408S-708N, with relatively weak temperature responses over the Southern Ocean and subpolar North Atlantic that even incudes cooling in places.nDIC is reduced under warming over much of the low latitudes but on average is larger over much of the Southern Ocean and also over the subpolar North Atlantic and Arctic.There is a strong correspondence between the patterns of SST and nDIC perturbations, with the area-weighted Spearman rank coefficient for the annual mean of both fields is 0.74.Globally, the average change in sea surface nDIC concentrations is 25.5% (28.3 mmol kg 21 8C 21 ).The average change over 458S-458N is 28.0%, over 908S-458S it is 11.8%, and over 458-908N it is 17.5%.Over 458S-458N the average change in surface nDIC is very similar to the net reduction in cumulative uptake over 458S-458N through air-sea fluxes (28.3%).The evolution of CO 2 fluxes integrated over 458S-458N in Fig. 2c is consistent with the interpretation that SST-driven decreases in the solubility of CO 2 of surface waters [K 0 5 CO 2 (aq)/pCO 2 ] as well as through perturbations to the dissociation constants for carbonic acid K 1 and K 2 (and thereby reflected in decreased nDIC) sustain decreased uptake of CO 2 for the warming case (COU) relative to the case without warming (BGC).
Positive changes in pCO 2 between COU and BGC (Fig. 3c) indicate a reduced thermodynamic gradient of pCO 2 into the ocean from the atmosphere.Aside from the Arctic, large positive pCO 2 changes are found in the subpolar North Atlantic and the Weddell Sea (Fig. 3c), with these changes reflected in globally maximum decreases in CO 2 uptake for the same regions (Fig. 3d).However, over most of 458S-458N, there is a smaller but widespread reduction in CO 2 uptake for COU relative to BGC, with exceptions being the equatorial Pacific upwelling regions and the regions spanning the Agulhas, the South Atlantic, and the Arctic.
In contrast to the spatially smoothly varying perturbation patterns for nDIC over the latitudes 458S-458N (Fig. 3b), the perturbation patterns for pCO 2 (Fig. 3c) and CO 2 fluxes (Fig. 3d) for COU relative to BGC reveal significantly more spatial structure.In fact, the eastern equatorial Pacific upwelling region exhibits locally enhanced uptake via gas exchange for COU relative to BGC, despite the fact that nDIC is reduced for this same upwelling region.For the case of the equatorial Pacific, this reflects the signature of reemergence of higher anthropogenic carbon concentrations for BGC relative to COU (Toyama et al. 2017;Zhai et al. 2017), with the meridional width of this extremum near the equator related to the air-sea equilibration time scale.
We consider in Fig. 3e the perturbation to the surface ocean Revelle factor (Revelle and Suess 1957) for the difference between COU and BGC averaged over 2080-99.Clearly the Revelle factor perturbations over large scales largely reflect the nDIC perturbations (Fig. 3b), with this in turn reflecting the SST perturbations (Fig. 3a).To emphasize the insight provided through nDIC, we also consider the change in DIC itself (Fig. 3f), where the differences between COU and BGC show significantly less correlation to the pattern of SST perturbations (Fig. 3a) than nDIC (Fig. 3b) due to the signature of freshwater perturbations under climate change on surface DIC concentrations.
What is particularly instructive about the differences shown in Fig. 3 is that COU exhibits higher sea surface pCO 2 and absorbs less CO 2 from the atmosphere over 458S-45, despite COU having lower nDIC and DIC concentrations than BGC.Although this may appear at first to be counterintuitive, we wish to emphasize here that it is the much higher DIC and thus pCO 2 for COU relative to PIN that enhances the sensitivity of pCO 2 to temperature changes, and thereby sustains the differences in CO 2 uptake over 458S-458N between COU and BGC.Likewise, it may seem counterintuitive that the Revelle factor is lower for COU than for BGC (Fig. 3e), while CO 2 uptake is less than COU relative to BGC (Fig. 3d).There as well, it is the impact of the much larger changes in DIC for COU relative to PIN that sustain the heat-carbon feedback mechanism.Whereas the perturbation to the Revelle factor averaged over 458S-458N between COU and BGC (Fig. 3e) is 20.6, the Revelle factor perturbation over the same region between COU and PIN is 4.0, with this being more than a factor of 6 larger in terms of absolute amplitude.
Our analysis in Fig. 3 focused on decadal mean perturbations for the model.Next, we turn our attention to seasonal variations in air-sea CO 2 fluxes and pCO 2 .Our interest is in the subtropical regions, and we begin in Fig. 4a with the example of monthly mean CO 2 fluxes integrated over 158-458N for the COU (black lines) and in Fig. 4b with monthly mean CO 2 fluxes integrated over the same region for BGC.For both cases, the seasonal cycle in CO 2 fluxes increases with time (consistent with Rodgers et al. 2008), with this increase being larger for COU than for BGC.To focus on the structure of the seasonal cycle at the end of the twenty-first century, we show next climatological fluxes over 2080-99 for 158-458N in Fig. 4c (COU shown in black; BGC shown in red).For this latitude range, there is a seasonal asymmetry in the amplification of amplitude for COU relative to BGC, with the amplitude perturbation being larger in late summer than during other parts of the seasonal cycle.Viewing BGC as a neutral state against which COU represents a perturbation, this clearly reveals that perturbations that are not evenly distributed over the seasonal cycle are sustaining the decreased uptake of COU relative to BGC of ;0.2 PgC yr 21 averaged over 2080-99.
For the latitude band 458-158S (Fig. 4d) the case of COU also reveals perturbations to the seasonal cycle of the neutral state BGC, with COU having ;0.1 PgC yr 21 less uptake over this latitude band averaged over the period 2080-99 than BGC.Once again, the difference between COU and BGC is not evenly distributed over the seasonal cycle.To get a foothold on the spatial structures over which this occurs, we consider for the North Pacific (the subtropical region with the largest signal) the late-summer DpCO 2 (ocean minus atmospheric pCO 2 ) averaged over August and September for COU (Fig. 4e) and BGC (Fig. 4f), both considered for climatologies constructed over 2080-99.Diagnostics of the piston velocities used for CO 2 gas exchange reveal that for COU they are only 2% smaller than for BGC over 458S-458N, implying that CO 2 uptake differences are dominated by surface ocean pCO 2 differences.The DpCO 2 patterns reveal that the thermodynamic gradient sustaining outgassing is significantly larger in its maximum in the subtropics and also characterized by a broader meridional structure for COU relative to BGC.
Next we consider the simulated seasonal evolution during the 2090s of pCO 2 within a phase space of SST and nDIC [analogous to that was presented in Fig. A1 of Nakano et al. (2015) and Fig. 5a of Schlunegger et al. (2019)] for the following four regions: 1) the formation region for NPSTMW (1408-1608E, 308-358N), 2) EQDIV (1408-1308W, 58S-58N), 3) the subpolar North Atlantic (SPNA) (408-308W, 508-608N), and 4) the Weddell Sea region of the Southern Ocean (WEDD) (408-308W, 628-728S).These regions are denoted with the bold black box in the four panels of Fig. 5. Representing the full seasonal evolution of pCO 2 allows us to assess the degree to which carbon-climate feedbacks are active in winter versus summer months, or whether they find expression over the full seasonal cycle.The first two of these (NPSTMW and EQDIV) lie within the latitudes spanned by the subtropical cells (458S-458N), whereas the other two regions (WEDD and SPNA) are important in the model for CO 2 exchange with the deeper ocean.(A comparison for the same locations against observational constraints for the 2000s is shown in Fig. 8.) In the NPSTMW region, each of the four simulations (COU, RAD, BGC, and PIN) presents a similar evolution of the seasonal cycle in pCO 2 in its absolute variations for SST and nDIC (Fig. 6a).COU_CC is not shown for the sake of clarity, as it closely resembles the case of COU.The annual mean nDIC offset between BGC and PIN is 238 mmol kg 21 , and between BGC and COU it is 13 mmol kg 21 .This decrease in nDIC for COU relative to BGC occurs within the context of an annual mean increase of SST of 2.38C, and a decrease in salinity and salinity-normalized total alkalinity of 0.25 psu and 2 mmol kg 21 , respectively, for COU relative to BGC.Importantly, the effect of salinity and nALK changes on pCO 2 is very small, and the offset between COU and BGC in SST-nDIC phase space occurs largely parallel to isolines of constant pCO 2 , albeit with the annual mean pCO 2 for COU being 11 matm higher than that of BGC, with this being only 2% of the total change in pCO 2 between BGC and PIN of 567 matm.We have seen with the case of COU_CC (Fig. 2c) that over the large scales of the subtropical cells the surface carbon cycle response in COU relative to BGC is largely due to SST perturbations.Viewed locally for STMW, Fig. 6a reveals that the amplitude of the nDIC response is approximately 6 mmol kg 21 8C 21 , with this being the approximate slope of the constant pCO 2 isolines shown in Fig. 6a.It is important to emphasize that although nDIC decreases by 5.3%, pCO 2 increases by 2.0%.Importantly, this response in a formation region for interior thermocline waters is smaller than the ;8% increase one would expect from SST-driven increases in pCO 2 in surface waters (the 4% pCO 2 increase per degree Celsius warming described by Weiss et al. 1982).
Next, we consider the case of EQDIV (Fig. 6b) for the 2090s, corresponding to the narrow strip of upwelling with 58 latitude of the equator between 1408 and 1308W.We have chosen to consider the EQDIV region as it offers a perturbation outgassing signature, thus being anomalous from the perturbation pCO 2 and CO 2 fluxes, while nevertheless being consistent with the broader nDIC perturbations occurring over the Pacific subtropical cells.We present analysis for a box with relatively narrow zonal extent to avoid the substantial zonal gradients in the broader Niño-3 box (1508-908W, 58S-58N), although the results considered here apply more generally for the region.Importantly, by focusing on a relatively narrow region, we target the upwelling zone for a midthermocline density class from the equatorial undercurrent.
The changes in annual mean pCO 2 between COU and BGC are very small (25 matm for EQDIV, as the annual mean local changes in nDIC (219 mmol kg 21 ) nearly compensate the annual mean changes in temperature (12.78C), while changes in salinity (20.05 psu) and nALK (21 mmol kg 21 ) have very small impacts.This near compensation in nDIC is reflected in a pCO 2 perturbation that is only 20.9% (a CO 2 ingassing perturbation) as large as the perturbation in pCO 2 for BGC relative to PIN for this same region in the 2090s, with the COU 2 BGC pCO 2 perturbation of 20.9% also being an order of magnitude smaller than what one would expect from thermodynamic temperature effects on pCO 2 .The equatorial Pacific upwelling region maintains a weak equatorial Pacific outgassing signature for both COU and BGC in the 2090s, but modestly weaker for COU.We interpret the lower nDIC concentrations for COU relative to BGC within the EQDIV region to not reflect reduced local uptake of CO 2 from the atmosphere, but rather to reflect differences in thermocline nDIC of recently upwelled waters along the equator.This is deduced from the existence of air-sea equilibration time scales in this region that exceed the residence time of surface  Weiss (1974) and K 1 and K 2 from Lueker et al. (2000).
waters in the EQDIV region.We interpret the reduced nDIC for COU relative to BGC in EQDIV to reflect rather differences in preformed nDIC in the extratropical source regions where thermocline waters form at the surface.
For the case of SPNA (Fig. 6c), nDIC is substantially higher for COU than it is for BGC.Consequently, pCO 2 is considerably higher (45 matm) for the COU case than for the BGC case over the full seasonal cycle, despite the compensating changes in SST (23.28C), salinity (21.91 psu), and nALK (117 mmol kg 21 ).These changes can be attributed to large-scale differences in the physical state of the ocean, and the response in biological processes in the ocean.As such, these changes are intrinsically linked to perturbations of the natural carbon cycle.This is followed by an analysis of the seasonal evolution in SST-nDIC space for pCO 2 for the case of WEDD (Fig. 6d).For COU considered as a perturbation to BGC, both nDIC and pCO 2 are higher (13 mmol kg 21 and 44 matm, respectively) with the increase in pCO 2 constituting a positive carbon-climate feedback.SST is cooler (20.58C) for COU relative to BGC, with this being most pronounced in summer, reflecting reduced convection in this region under global warming (Manabe et al. 1991).From this analysis we see that temperature perturbations that compensate for nDIC as well as salinity (20.38 psu) and nALK (11 mmol kg 21 ) are relatively minor.The perturbations to nDIC that drive the pCO 2 perturbations are largest in late austral winter (September).
Taken together, our analysis emphasizes the contrast between two regions within 458S-458N in the Pacific (Figs. 6a,b) and the two high-latitude regions (Figs.6c,d).The high-latitude regions exhibit strong perturbations in the natural carbon cycle, associated with changes in circulation and/or biology, consistent with behavior seen in previous modeling studies (Sarmiento et al. 1998;Winton et al. 2013;de Lavergne et al. 2014).For both highlatitude sites, nDIC increases for COU relative to BGC despite decreasing in temperature, implicating mixing and/or biological changes that perturb nDIC.The two sites within 458S-458N, on the other hand, are characterized by surface nDIC perturbations that are to first order driven by temperature perturbations, with gas exchange playing an important supporting role.The relatively structured pCO 2 perturbation field shown in Fig. 3c represents the interplay of spatial variations of the slope of pCO 2 lines in SST-nDIC space for the different regions of the surface ocean between COU and BGC.
We have identified broadscale reductions in response to surface warming of both sea surface nDIC and cumulative CO 2 uptake over 458S-458N (both by ;8%) by the end of the twenty-first century.The relatively smooth distribution of nDIC differences (Fig. 3b) over the equatorial upwelling region indicates that reemerged equatorial thermocline waters have nDIC and thermal perturbations that reflect those of the overlying surface waters.This is consistent with thermocline waters more generally having nDIC concentrations that are reduced by 8% due to heating perturbations to surface waters, in particular for the Eastern Subtropical Mode Water (ESTMW) formation regions thought to be important as a source of Equatorial Undercurrent (EUC) waters (Rodgers et al. 2003;Goodman et al. 2005).The fact that pCO 2 displacements for the EQDIV region in Fig. 6b between BGC and COU fall nearly along a constant pCO 2 isoline reflects large-scale shallow overturning processes (Nakano et al. 2015;Toyama et al. 2017;Zhai et al. 2017) rather than simply local processes in the equatorial divergence region.These results are consistent with the observationally based analysis of Fine et al. (2001) that revealed a decadal-scale shallow overturning time scale for thermocline waters.
By establishing over the region 458S-458N a specific mechanistic coupling between temperature and carbon, our results provide a mechanistic underpinning of the g oc factor (g oc being the carbon-climate feedback diagnostic of Friedlingstein et al. (2003), for the reduction in cumulative CO 2 uptake per degree Celsius warming).Importantly, the framework of invoking SST-nDIC phase space to interpret perturbations (Fig. 6), in conjunction with our earlier analysis in Figs. 3 and 4, has allowed us to identify that the feedback over 458S-458N is triggered by surface heating, and subsequently the sensitivity of surface pCO 2 to temperature perturbations is amplified in time through the invasion flux of anthropogenic carbon (taken as COU minus PIN), in particular during the warm season.Importantly, Fig. 6 also confirms that biological perturbations are not important players in sustaining the differences in CO 2 uptake between COU and BGC over 458S-458N.
The details of the fields used to identify g oc by Friedlingstein et al. (2003) are not equivalent to the local relationship we identified in Fig. 6 between SST and nDIC.Thus we are motivated to consider over the expanse 458S-458N the relationship between the SST perturbations (COU minus BGC) over this region with the perturbation to cumulative CO 2 uptake over this region (COU minus BGC) in Fig. 7a (analogous to Fig. 2d in Friedlingstein et al. 2006).For this analysis 10-yr means of both fields over 1861-2100 were used.As temperatures over 458S-458N increase, there is a corresponding decrease in CO 2 uptake for COU relative to BGC.The response of cumulative CO 2 uptake to warming is relatively weak until approximately the year 2000, with an abrupt change to a stronger sensitivity over the twenty-first Unauthenticated | Downloaded 04/08/24 11:43 AM UTC century.The transition in this sensitivity occurs at approximately 0.88C.In other words, there is a pronounced nonlinearity in the relationship between large-scale temperature perturbations and large-scale cumulative carbon uptake perturbations.Next cumulative CO 2 uptake perturbations over 458S-458N are considered against globally averaged surface air temperature (Fig. 7b) shows a similar structure, reflecting the close correspondence of SSTs averaged over 458S-458N with globally averaged surface air temperature (Fig. 7c).Taken together, this analysis suggests that the temporal structure of cumulative ocean CO 2 flux perturbations to global surface temperature perturbations presented by Friedlingstein et al. (2006) are reflecting the heat-carbon feedback mechanism with an underlying nonlinear response.Last but not least, the relationship between perturbations to cumulative CO 2 uptake over 458S-458N are shown against changes in the mean nDIC concentrations (COU minus BGC) averaged over 458S-458N (Fig. 7d).This reveals a close correspondence between two fields, but with a similar curvature to what is seen in Fig. 7a.
As a complement to the analysis presented in Fig. 6 for the projected 2090s evolution of pCO 2 in SST-nDIC phase space, we consider in Fig. 8 a similar analysis for the modeled behavior of pCO 2 over the 2000s (specifically a climatology over 2000-09) against two climatologies derived from observations over the same period.In particular, we consider variability over three of the four regions that were presented in Fig. 5, using the interannually varying monthly gridded products of Majkut et al. (2014) (PU-MCMC) and Iida et al. (2015) (JMA).For three of these regions, the modeled evolution FIG. 7. (a) Relationship between perturbations to SST perturbations (8C) averaged over 458S-458N (COU minus BGC) taking 10-yr means between 1861 and 2100 and perturbations to cumulative CO 2 uptake (PgC) (COU minus BGC) over 458S-458N; (b) relationship between perturbations to globally averaged surface air temperature (8C) (COU minus BGC) (8C) and perturbations to cumulative CO 2 uptake (PgC) (COU minus BGC) over 458S-458N; (c) relationship between perturbations to global atmospheric surface air temperature (TAS) perturbations (8C) averaged over 458S-458N (COU minus BGC) taking 10-yr means between 1861 and 2100 and perturbations to SST (8C) averaged over 458S-458N (COU minus BGC) taking 10-yr means between 1861 and 2100; (d) relationship between perturbations to cumulative CO2 uptake (PgC) over 458S-458N and perturbation to spatially integrated surface nDIC (mmol kg 21 ) over 458S-458N.
for COU over the 2000s (ESM2M_COU) is evaluated against the observational constraints.Although the ESM2M_ COU evolution is also shown for the WEDD region in the 2000s (Fig. 8d), data constraints there were deemed to be insufficient to merit comparison against the modeled state.
For the three cases where the model is explicitly compared with the observations (Figs.8a-c), the two observational products (PU_MCMC and JMA) are in fact more similar to each other than to the model (ESM2M_COU), supporting our assumption that the observational products are valid for assessing model biases.For the NPSTMW site (Fig. 8a) the analysis indicates that agreement between all three estimates is best in winter, with model biases of order 10 mmol kg 21 for nDIC and of order 18C for SST, and biases toward a weak pCO 2 of order 10-20 matm.Biases in summer are approximately twice as large for this region.For the case of EQDIV (Fig. 8b) the model exhibits a bias with nDIC that is low by 10-20 mmol kg 21 , and an annual mean pCO 2 that is too low by 10-40 matm, although the amplitude in the seasonal cycle in nDIC and SST well approximates those identified in the observational products.For the subpolar North Atlantic (Fig. 8c) the model reveals a substantial bias (too warm by 28-38C and nDIC concentrations too low by 10-30 mmol kg 21 , and pCO 2 in winter (the time of minimum temperatures) is 25-50 matm smaller than the observational products.For the Weddell Sea, the evolution of pCO 2 in SST-nDIC phase space is shown for the model for consistency, but as was stated earlier data sparsity issues in these regions render these fields to be inappropriately represented in the PU_ MCMC and JMA observationally based products.

Discussion
We set out to investigate the degree to which heat and carbon coupling under enhanced anthropogenic carbon emissions can occur through the impact of heat on the For WEDD data constraints were deemed to not be sufficient to be included, but the model behavior is shown.The figure uses K 0 from Weiss (1974) and K 1 and K 2 from Lueker et al. (2000).solubility of CO 2 and ensuing perturbations to the CO 2 buffering capacity of seawater under CO 2 invasion, in the absence of perturbations to the physical circulation state or the biological state of the ocean.We began with a global analysis, and we were able to identify that nearly half (48%) of the global ocean carbon-climate feedback occurs over the latitude range 458S-458N for GFDL's ESM2M, with approximately 60% of the global uptake of anthropogenic CO 2 occurring over this same latitude range.To attribute the degree to which this behavior over 458S-458N is due to solubility effects rather than perturbations to the ventilation or biological state of the ocean, we devised a new simulation (COU_ CC) that maintained a preindustrial physical circulation state of the ocean, but includes the impact of climate warming from our COU simulation.This was accomplished through the representation of temperature and salinity in the calculation of pCO 2 exclusively in the gas exchange routines of the model.With this COU_CC simulation, we were able to confirm that it is the impact of warming on solubility, that is enhanced by the diminished buffering capacity of CO 2 in seawater in a high-CO 2 world, rather than changes in ventilation or stratification, that dominate the marine carbon-climate feedbacks over 458S-458N.
The analysis of cumulative carbon perturbations (COU minus BGC) versus SST perturbations (COU minus BGC) over 458S-458N in Fig. 7a illustrated a sensitivity that increases with time, in particular above the 18C warming perturbation threshold.This delayed response in cumulative CO 2 uptake to warming may well reflect the decadal water mass renewal time scales for the shallow overturning, in particular an intergyreexchange time scale for thermocline waters that upwell along the equator (Fine et al. 2001).In other words, as the offset of COU from BGC for EQDIV (Fig. 6b) is not a local but rather a delayed response to conditions in the extratropical source regions (upwelling of waters that have had their CO 2 buffering capacity reduced through the invasion flux of CO 2 from the atmosphere) related to the equatorial region through thermocline transport.
To reiterate our principal result, the heat-carbon feedback is triggered over 458S-458N by surface warming, and subsequently the sensitivity of sea surface pCO 2 to further temperature perturbations is amplified by the cumulative invasion flux of C ant and its impact on surface DIC concentrations.It is then important to consider the degree to which the results presented here may be model-specific.One area of concern with ESM2M is that it has been documented to be on the low end of the CMIP5 models in terms of its transient climate sensitivity, meaning that feedbacks could be weaker with this model than other CMIP5-class models.Even with this model expected to experience a weaker TCR (Forster et al. 2013), including a weaker warming response over the low latitudes, the heat-carbon feedback identified here is nevertheless interpreted to be robust.The question also arises as to whether the model's temperature bias in summer in Fig. 8a (the North Pacific STMW formation region) complicates our interpretation of the importance of the warm season in sustaining the heatcarbon feedback.There we wish to emphasize that the amplitude of the seasonal cycle in SST in the model has a bias of less than 10%, which is not large.It is in fact the model skill with the seasonal amplitude that is of particular interest for our inferred mechanism.It also warrants mention that the mapping uncertainties for carbon-related variables within observational products are larger than 10%, as is clear from a number of studies (Wanninkhof et al. 2013;Rödenbeck et al. 2015;Gregor et al. 2019), with biases also being large for the seasonal cycle on local scales.We recommend that future work address seasonal biases seen here in ESM2M comprehensively more broadly in CMIP-class models.
Recently there has been interest in evaluating Earth system models to develop emergent constraints to identify clues from historical records to constrain projected future changes in the Earth system (e.g., Cox et al. 2013).It warrants mention here that there is a stark contrast for the equatorial Pacific region between the climatological seasonal changes emphasized in Fig. 6b and the variations of pCO 2 on El Niño time scales for the model (not shown).The amplitude of the climatological seasonal cycle in pCO 2 is small, with displacements in SST-nDIC space largely along the derived stationary pCO 2 isolines, whereas for El Niño-Southern Oscillation (ENSO) time scales the variability is large with important displacements of pCO 2 in SST-nDIC space that are transverse to the derived stationary pCO 2 isolines in the figure.In the language of Takahashi et al. (2002) and Feely et al. (2002), the feedback as manifested in pCO 2 over the EQDIV region (the offset of COU relative to BGC) is SST driven, with this being confirmed over larger scales with the COU_CC simulation shown in Fig. 2.This then stands in contrast to ENSO time scales, where the changes in pCO 2 in the equatorial Pacific upwelling region were identified by Takahashi et al. (2002) and Feely et al. (2002) as being DIC driven (circulation and biology) rather than SST driven.We interpret this to provide a note of caution regarding the consideration of ''emergent constraints'' as a means to relate observed changes in the modern ocean state to project future changes.Although this has been argued to work for marine net primary productivity (NPP) (Kwiatkowski et al. 2017) and for terrestrial carbon-climate feedbacks (Cox et al. 2013), the mechanisms controlling carbon feedbacks are distinct here from the mechanisms controlling ENSO variability in sea surface pCO 2 .
Although our principal focus in this study has been on the lower latitudes, there are two points with respect to the higher latitudes that warrant mention.First, as we pointed out for the subpolar North Atlantic (Fig. 6c) and the Weddell Sea (Fig. 6d), the offset of the COU and BGC runs relative to the constant pCO 2 isolines in the figures indicates changes that are primarily DIC driven rather than SST driven.In fact, this holds over larger scales, as is revealed for the integrated CO 2 fluxes over the high latitudes (Fig. 2b).There the integrated CO 2 fluxes for COU_CC more closely follow BGC than COU.The second point that we wish to emphasize for the high latitudes is that the feedbacks there may become important after they have become important over the low latitudes.This can be inferred by comparing the degree of separation between COU and BGC over the high latitudes (Fig. 2b) and the low latitudes (Fig. 2c).Resolving this more quantitatively will require a large ensemble approach for both BGC and COU, as a means to distinguish between forced responses and natural variability.
The mechanism emphasized here over the low latitudes is consistent with the connection made between water column temperature perturbations under climate change and carbon-climate feedbacks identified by Schwinger and Tjiputra (2018), albeit with different perturbations applied to different Earth system models.Although well beyond the scope of this study, a constructive path forward may well be to extend the analyses with CMIP5 and CMIP6 models presented by Arora et al. (2020) for which the array of COU and BGC simulations exist with 1% per year atmospheric forcing for quantifying feedbacks by region considered here and in Table 2 of Roy et al. (2011).This will help to resolve the following question: Do CMIP5 models tend to exhibit a stronger low-latitude contribution to global carbon-climate feedbacks than ESM2M (the model considered here), with this reflecting the weak thermal transient climate response found for ESM2M?
The analyses of pCO 2 evolution in SST-nDIC phase space in Figs. 6 and 8 does not in itself reveal whether and how modulations of the seasonal cycle resulting from changes in the CO 2 buffering capacity of seawater may contribute to carbon-climate feedbacks.These results are consistent with the findings of the earlier studies (Rodgers et al. 2008;Riebesell et al. 2009;Gorgues et al. 2010;Hauck and Völker 2015;Landschützer et al. 2018;Fassbender et al. 2018) that have emphasized the growth of the seasonal cycle in sea surface pCO 2 under the joint effects of the invasion flux of C ant and warming.It is left as a subject for future investigation to explore the degree to which carbon-climate feedbacks over the low latitudes emphasized here are stronger over summer than winter months.

Conclusions
Our analysis has identified a heat-carbon feedback mechanism that sustains approximately half of the global carbon-climate feedback by 2100 under historical/RCP8.5 perturbations with GFDL's ESM2M Earth system model.Over the latitude range 458S-458N, where approximately 60% of the global ocean uptake of anthropogenic carbon (C ant ) occurs in ESM2M, the heat-carbon feedback emphasized here is dominant.Importantly, the heat-carbon feedback is distinct from and weaker than what one would expect from the effect of solubility (through surface ocean warming) on pCO 2 alone.Rather, this low-latitude feedback is initially triggered by surface temperature perturbations, but it is then further amplified through the effect of the net cumulative invasion flux of anthropogenic carbon into the surface ocean.This in turn enhances the sensitivity of pCO 2 temperature changes, in particular during summer conditions.As a consequence of less CO 2 uptake by the ocean in a higher CO 2 world, more CO 2 is left in the atmosphere, enhancing further warming.
A rectification of the seasonal cycle in pCO 2 in the subtropics was identified to play an important role in sustaining the heat-carbon feedback identified here.This occurs in particular through enhanced summer outgassing of CO 2 for COU relative to BGC in the subtropics, but also more generally through a modulation of the full seasonal cycle.The fact that seasonal pCO 2 variability in the subtropics is amplified for COU relative to BGC is due to the full invasion flux (COU minus PIN) of C ant , with this overwhelming the much smaller compensation effect seen in Revelle factor differences between COU and BGC.The fact that the feedback over 458S-458N is relatively weak for perturbation temperatures less than 0.88C, and then increases strongly for temperature perturbations greater than 0.88C (Fig. 7a) is interpreted to represent a delay associated with the renewal time scales for thermocline waters (Fine et al. 2001).This serves to underscore the importance of reemergence of C ant into the ocean's mixed layer as a mechanism that can modulate the strength of the heat-carbon feedback.
The introduction of the COU_CC run (invoking data override) provided independent support for our main mechanistic interpretation, namely that the heat-carbon feedback dominates carbon-climate feedbacks over 458S-458N.The experiment confirms that the feedback is triggered by heating and then amplified by the invasion flux of C ant , rather than being driven by ventilation rate perturbations (circulation state changes) or perturbations to biology over 458S-458N.This study complements previous works that have emphasized positive marine carbon-climate feedbacks that are driven by the combined effects of circulation state and biological perturbations over the high latitudes (Maier-Reimer et al. 1996;Sarmiento et al. 1998).It is our hope that this study will motivate future explorations of how the heat and carbon coupling mechanisms discussed here operate over higher latitudes, where there is a more complicated interplay between a broader variety of mechanisms to sustain feedbacks.

FIG. 3 .
FIG. 3. Maps showing annual mean perturbation structures, averaged over the 20-yr intervals 2080-99.The fields shown are (a) difference in sea surface temperature (8C) between COU and BGC; (b) difference in sea surface nDIC (mmol kg 21 ) for COU and BGC runs; (c) difference in pCO 2 (matm) between COU and BGC; (d) difference in sea-air CO 2 flux (molC m 22 yr 21 ) between COU and BGC; (e) the change in the surface ocean Revelle factor between COU and BGC; (f) difference in sea surface DIC (mmol kg 21 ) for COU and BGC runs.