Amplified Subsurface Signals of Ocean Acidification

We evaluate the impact of anthropogenic carbon (Cant) accumulation on multiple ocean acidification (OA) metrics throughout the water column and across the major ocean basins using the GLODAPv2.2016b mapped product. OA is largely considered a surface‐intensified process caused by the air‐to‐sea transfer of Cant; however, we find that the partial pressure of carbon dioxide gas (pCO2), Revelle sensitivity Factor (RF), and hydrogen ion concentration ([H+]) exhibit their largest responses to Cant well below the surface (>100 m). This is because subsurface seawater is usually less well‐buffered than surface seawater due to the accumulation of natural carbon from organic matter remineralization. pH and aragonite saturation state (ΩAr) do not exhibit spatially coherent amplified subsurface responses to Cant accumulation in the GLODAPv2.2016b mapped product, though nonlinear characteristics of the carbonate system work to amplify subsurface changes in each OA metric evaluated except ΩAr. Regional variability in the vertical gradients of natural and anthropogenic carbon create regional hot spots of subsurface intensified OA metric changes, with implications for vertical shifts in biologically relevant chemical thresholds. Cant accumulation has resulted in subsurface pCO2, RF, and [H+] changes that significantly exceed their respective surface change magnitudes, sometimes by >100%, throughout large expanses of the ocean. Such interior ocean pCO2 changes are outpacing the atmospheric pCO2 change that drives OA itself. Re‐emergence of these waters at the sea surface could lead to elevated CO2 evasion rates and reduced ocean carbon storage efficiency in high‐latitude regions where waters do not have time to fully equilibrate with the atmosphere before subduction.


Introduction
Anthropogenic carbon (C ant ) in the atmosphere is primarily transferred to the ocean through the air-sea exchange of carbon dioxide gas (CO 2 ).To date, the ocean has absorbed ∼25% of the total anthropogenic emissions since industrialization (Friedlingstein et al., 2022).Persistent growth in the atmospheric partial pressure of CO 2 (pCO 2 ) has caused persistent growth in upper ocean C ant concentrations, with the highest values generally occurring at the sea surface.However, the accumulation of C ant is not spatially uniform (Gruber et al., 2019;Sabine et al., 2004).Greater penetration of C ant in deep water formation regions and the subtropics is caused by regional circulation that transports surface waters to the ocean interior (DeVries, 2014;Iudicone et al., 2016;Khatiwala et al., 2013).
In addition to circulation, heterogeneity in surface ocean carbonate chemistry influences the impact of increasing atmospheric pCO 2 on dissolved inorganic carbon (DIC), and thus the efficiency of ocean C ant accumulation (Revelle & Suess, 1957).Meridional surface gradients in temperature and salinity, in addition to physical and biological processes that can modulate local air-sea CO 2 disequilibrium, dictate the meridional surface gradient in DIC concentration (Hamme et al., 2019;Volk & Hoffert, 1985) and the associated DIC-to-total alkalinity ratio (Wu et al., 2019); a measure of the carbonate system buffer capacity (Egleston et al., 2010).The carbonate system Abstract We evaluate the impact of anthropogenic carbon (C ant ) accumulation on multiple ocean acidification (OA) metrics throughout the water column and across the major ocean basins using the GLODAPv2.2016bmapped product.OA is largely considered a surface-intensified process caused by the air-to-sea transfer of C ant ; however, we find that the partial pressure of carbon dioxide gas (pCO 2 ), Revelle sensitivity Factor (RF), and hydrogen ion concentration ([H + ]) exhibit their largest responses to C ant well below the surface (>100 m).This is because subsurface seawater is usually less well-buffered than surface seawater due to the accumulation of natural carbon from organic matter remineralization.pH and aragonite saturation state (Ω Ar ) do not exhibit spatially coherent amplified subsurface responses to C ant accumulation in the GLODAPv2.2016bmapped product, though nonlinear characteristics of the carbonate system work to amplify subsurface changes in each OA metric evaluated except Ω Ar .Regional variability in the vertical gradients of natural and anthropogenic carbon create regional hot spots of subsurface intensified OA metric changes, with implications for vertical shifts in biologically relevant chemical thresholds.C ant accumulation has resulted in subsurface pCO 2 , RF, and [H + ] changes that significantly exceed their respective surface change magnitudes, sometimes by >100%, throughout large expanses of the ocean.Such interior ocean pCO 2 changes are outpacing the atmospheric pCO 2 change that drives OA itself.Re-emergence of these waters at the sea surface could lead to elevated CO 2 evasion rates and reduced ocean carbon storage efficiency in high-latitude regions where waters do not have time to fully equilibrate with the atmosphere before subduction.

Plain Language Summary
The chemistry of the upper ocean is changing due to the absorption of excess carbon in the atmosphere resulting from human activities, a process commonly referred to as ocean acidification (OA).The highest concentrations of excess carbon in the ocean are generally found at the surface; however, some of the largest chemical changes resulting from this carbon buildup are occurring below the sea surface.This subsurface intensification is caused by nonlinear responses of some chemical parameters to opposing vertical gradients of natural carbon versus excess carbon.Our findings emphasize the need to study multiple metrics of OA throughout the water column to comprehensively evaluate changes in the habitability of interior ocean realms and potential connections to ocean carbon storage timescales.
FASSBENDER ET AL.For a given relative change in pCO 2 , a low RF indicates that a relatively large change in DIC will occur, and a high RF indicates that a relatively small change in DIC will occur.The mean surface ocean RF is ∼10, indicating that a 10% change in pCO 2 is required to elicit a 1% change in DIC.However, spatial variability in the RF indicates that some areas of the ocean are more efficient than others at accumulating DIC for a given surface ocean pCO 2 perturbation (Figure 1a and Figure S1 in Supporting Information S1; Fassbender et al., 2017).Below the surface, the RF increases with depth, more than doubling in many ocean regions, indicating an increase in the relative sensitivity of pCO 2 to DIC (Figure 1c).
In the ocean interior, the three-dimensional habitats of marine species are shaped by environmental conditions (e.g., Deutsch et al., 2020;Howard et al., 2020) that evolve with natural variability and external forcing.Yet, few studies have evaluated how OA metrics influence the distributions of economically important species and their prey (e.g., Guinotte & Fabry, 2008).For example, Ω Ar is a valuable indicator of the energetic requirement to sustain net calcification (e.g., Bach, 2015) for species that may serve as prey for pelagic fish or be economically important in their own right (Cooley & Doney, 2009).Additionally, there is evidence that some fish species may be severely impacted by elevated seawater pCO 2 levels that can make it difficult to expel CO 2 from the gills, causing a buildup of CO 2 in the blood (i.e., hypercapnia; Heuer & Grosell, 2014).Understanding how OA is evolving below the surface, where new mesopelagic fisheries may become important to feed the growing global population (Fjeld et al., 2023;Kourantidou & Jin, 2022), is a critical first step toward identifying potential downstream impacts of OA on marine species distributions and ecosystem health.
Building on prior studies, we evaluate the impact of C ant accumulation on multiple commonly considered OA metrics throughout the water column and across the major ocean basins.We then quantify the impact of carbonate system nonlinearities on the C ant -induced OA metric changes.The indirect impacts of C ant on ocean circulation and the biological carbon pump through climate change are additional complications that we do not address in this study.Finally, we consider potential implications for interior ocean habitability, carbon storage efficiency, and impacts resulting from marine carbon dioxide removal strategies.10.1029/2023GB007843

Methodology
3 of 14 ard depth levels and normalized to the year 2002 as an annual climatology: the GLODAPv2.2016bmapped product (Key et al., 2015;Lauvset et al., 2016).Salinity (S), temperature (T), total alkalinity (TA), DIC, and nutrient (N; phosphate and silicate) data were used to compute pH, Ω Ar , pCO 2 , [H + ], and RF (Figure S2 in Supporting Information S1).All carbonate system calculations in our study were made in MATLAB using the CO2SYSv3 (version 3.2.0;Sharp et al., 2020) carbonate system calculator with the following user settings: total pH scale; equilibrium constants of Lueker et al. (2000); sulfate dissociation constant of Dickson et al. (1990); fluoride dissociation constant of Perez and Fraga (1987); and borate-to-salinity ratio of Lee et al. (2010).
C ant values included in the GLODAPv2.2016bmapped product were subtracted from the DIC values to estimate a quasi-preindustrial (PI) DIC concentration (DIC PI ).The GLODAPv2.2016bmapped product does provide a DIC PI field, but it is mapped from calculations of DIC PI at observation locations and is nearly identical (0.0 ± 0.2 μmol kg −1 ) to the mapped DIC minus mapped C ant .We chose to deal with DIC PI calculated from the mapped DIC and C ant for consistency.Salinity, temperature, TA, DIC PI , and nutrient data were then used to compute pH PI , Ω Ar PI , pCO 2 PI , [H + ] PI , and RF PI .This approach isolates the influence of accumulated C ant in setting ocean chemical conditions (in 2002 c.e.) but does not account for changes in temperature, salinity, ocean circulation, or biology that may have occurred since industrialization.Differences between the year 2002 and PI values are used to evaluate subsurface patterns of OA metric changes caused by accumulated C ant (for a variable X this would be  Δ   ; Figure S3 in Supporting Information S1), with an example calculation for pH shown here: We consider subsurface changes in each OA metric relative to the respective mean change within the upper 50 m of the local overlying water ( Local overlying waters do not typically reflect the origin conditions of the waters below; however, using the local, rather than global, mean change in near-surface waters as a reference helps to reveal how the influence of this feedback varies with depth (Figures S4 and S5 in Supporting Information S1).

Nonlinear Component of OA Metric Changes Caused by C ant
To quantify how carbonate system nonlinearities may be contributing to the total change in each OA metric, we use previously published, globally mapped estimates of preformed properties (Carter et al., 2021) provided on the same grid as the GLODAPv2.2016bmapped product.These so-called preformed properties reflect the property concentrations that would be expected in the absence of carbon accumulated through organic matter remineralization (Figure S6 in Supporting Information S1) and calcium carbonate dissolution (Figure S7 in Supporting Information S1).The interior ocean preformed property estimates rely on the transport matrix output from a data-assimilation ocean circulation inverse model (OCIM: DeVries, 2014) to determine where interior ocean water masses were last near the ocean surface, and two alternative circulation products are used as a component of the uncertainty estimate for the product.OCIM in turn relies on GLODAP observations (Key et al., 2004) of potential temperature, salinity, radiocarbon, and CFC-11 to estimate the climatological mean state ocean circulation.Locally interpolated regression algorithms developed by Carter et al. (2017), which are also based on GLODAP observations, were used to estimate preformed property values at the base of the mixed layer.These values were propagated into the ocean interior using the OCIM transport matrix to estimate the interior ocean preformed property values.As a final step, Carter et al. (2021) regridded the preformed properties to match the GLODAPv2.2016bmapped product grid for public release of the files.
We focus on the biogenic component of natural carbon in the ocean due to its dominance in setting the vertical carbon gradient (DeVries, 2022;Sarmiento & Gruber, 2006) and because it is the addition of biogenic carbon, after natural and anthropogenic carbon have been added through air-sea exchange, that induces the nonlinear subsurface changes we aim to characterize (Figure S8 in Supporting Information S1).We subtract the preformed dissolved oxygen (O 2 0 ) field from the GLODAPv2.2016bmapped product O 2 field to estimate the true oxygen utilization (TOU).To convert TOU to respiratory DIC, we use the canonical respiratory quotient (O 2 :C) of (−170/117) that was estimated from mixed phytoplankton tows (Hedges et al., 2002) and interior ocean water mass analyses (Anderson & Sarmiento, 1994) to evaluate the bulk remineralized carbon pool in the ocean interior (sensitivity analysis presented in Text S1 and Figures S9-S10 in Supporting Information S1).We subtract preformed nitrate, silicate, and phosphate from the GLODAPv2.2016bmapped product fields for each parameter to determine their respiratory burdens.Following Carter et al. (2021), we compute potential alkalinity (pTA) by subtracting the GLODAPv2.2016bmapped product nitrate field multiplied by 1.36 (Wolf-Gladrow et al., 2007) from the TA field to account for the minor influence of remineralization on TA.We perform the same calculation using the preformed TA (TA 0 ) and preformed nitrate fields to estimate preformed potential TA (pTA 0 ).The difference between pTA and pTA 0 yields the calcium carbonate dissolution term (TA CaCO3 ).Subtracting TA 0 and TA CaCO3 from the TA field yields the residual, remineralized TA component.We multiply TA CaCO3 by 0.5 to calculate the DIC calcium carbonate dissolution component.Preformed DIC (DIC 0 ) is then calculated by subtracting the respiratory and calcium carbonate DIC components from the GLODAPv2.2016bmapped product DIC field.
For each OA metric, the preindustrial preformed (PI,0) value is subtracted from the year 2002 preformed value to determine OA metric changes  ( Δ ) in the absence of natural biogenic byproducts, with an example calculation for pH shown here: For each OA metric (e.g., X), we then evaluate the difference between C ant -induced changes occurring with (  Δ  ant ) and without (  Δ  ant 0 ) the background of natural biogenic carbon byproducts (hereafter referred to as natural carbon).This allows us to quantify how much of the OA metric response is caused by nonlinear carbonate chemistry effects induced by natural and anthropogenic carbon pool interactions (Figure S8 in Supporting Information S1), with an example calculation for pH shown here: In the history of a water mass, C ant accumulates before and during water mass formation and subduction, whereas natural carbon accumulates after water mass subduction.Therefore, this nonlinear effect should be considered the impact of C ant upon the magnitude of change induced by natural carbon rather than the reverse.Nevertheless, hereafter we discuss the impact in the reverse sense because we are interested in how human emissions modulate the signals associated with natural environmental variability.10.1029/2023GB007843 5 of 14 Our calculations rely on numerous assumptions and data products, each of which contains a degree of uncertainty.However, in Text S2, we describe a validation exercise based on independent output from a global ocean biogeochemistry model that reproduces the significant patterns that we highlight in our findings (Figure S11 in Supporting Information S1).This confirms that our findings are not spurious products of our methodological choices but are robust signals that are consistent with modern understanding of ocean circulation, biogeochemistry, and anthropogenic forcing.While this exercise shows that models already represent subsurface amplification of OA metric changes through carbonate system nonlinearities, our study serves to characterize and to call attention to this combination of unintuitive feedbacks and to provide observational evidence for its existence and extent in the real ocean.

Estimation of Uncertainties
The GLODAPv2.2016bmapped product provides spatially resolved error fields for each parameter.These errors reflect mapping uncertainties and do not include measurement uncertainties or calculation errors (for pH and Ω Ar ), as mapping errors are thought to dominate (Lauvset et al., 2016;Olsen et al., 2016).We compute the [H + ] error directly from the provided pH and pH error fields (Text S3).The Carter et al. (2021) preformed a property mapped product also provides spatially resolved error fields for each parameter.We use the error fields in Monte Carlo simulations and standard error propagation procedures (i.e., summing in quadrature) to estimate the uncertainty in the calculated values presented herein.For pCO 2 , RF, and all PI values, we perform a 1,000 iteration Monte Carlo simulation of the CO2SYSv3 calculations while individually varying all input values around their provided errors in a Gaussian manner using a MATLAB pseudorandom number generator (randn) that produces a set of numbers with a mean equal to zero and standard deviation equal to one.We also calculate the change in each parameter from the PI period to the year 2002 within each simulation.The standard deviation across the 1,000 simulations for each parameter is used to estimate its uncertainty.Standard error propagation is used in subsequent calculations as needed.The computed uncertainties are used for stippling in various figures to convey statistical confidence.

Identifying OA Metrics With Subsurface Intensified Changes
pH and Ω Ar do not display coherent, amplified subsurface C ant -induced changes  ( Δ  ant ) relative to the local mean changes within the upper 50 m along the transects evaluated (Figure 2 and Figure S12 in Supporting Information S1).However, there are patchy regions of slightly elevated subsurface  ΔpH   that exceed the calculation uncertainties (Figures 2a-2c), and similar signals have been identified previously (Carter et al., 2019;Lauvset et al., 2020).pCO 2 and [H + ] display amplified subsurface C ant -induced changes, relative to the local mean changes within the upper 50 m, with the largest changes occurring at depths with moderate to low C ant concentrations (Figure 2 and Figure S12 in Supporting Information S1).These features are spatially coherent and found throughout significant portions of each ocean transect, reflecting a signal that is global in nature.The largest subsurface signals occur in the North Pacific (∼100-500 m) where the waters are old (Figure S13 in Supporting Information S1), the vertical pCO 2 and DIC gradients are steep (Figure S14 in Supporting Information S1), and carbonate system buffering is naturally weak (Figure 1c).In these areas, small changes in DIC have a large influence on pCO 2 due to the elevated background pCO 2 :DIC ratio (Figure S15a S15b in Supporting Information S1).Where DIC < TA, as it is in much of the modern ocean, the RF increases with DIC additions, as carbonate ions are consumed to neutralize the added CO 2 (Figure S15c in Supporting Information S1), thus diminishing the quantity of carbonate ions remaining to buffer the seawater against further changes.When DIC > TA, the carbonate ions have been largely consumed, leading to an approximately linear increase in pCO 2 with increasing DIC (Figures S15a and S15c in Supporting Information S1).In these regions, the RF declines with added DIC, primarily due to a constant ΔpCO 2 /ΔDIC and rapidly growing pCO 2 denominator in Equation 1.Like pCO 2 , subsurface  ΔRF   values throughout the Southern Ocean are larger than the uncertainties (Figures 2m-2o).

Nonlinear Interactions Between Natural and Anthropogenic Carbon Pools
Regional differences in vertical DIC and TA gradients caused by heterogeneous accumulation of biogenic and anthropogenic carbon (Figures S6 and S7 in Supporting Information S1) lead to regional differences in the nonlinear responses of OA metrics to natural and anthropogenic carbon pool interactions (Figure 3).Ocean C ant distributions Nonlinear interactions between natural and anthropogenic carbon pools are working to both dampen and amplify OA metric changes in different regions of the water column.For pH, nonlinearities amplify the decrease in pH throughout coherent regions of the subsurface ocean (Figure 3); however, this effect is small, accounting for <25% of the For Ω Ar , carbonate system nonlinearities dampen the decrease in Ω Ar throughout the water column, but the effect is small, accounting for <25% of the total Ω Ar change in most waters above the 10 μmol kg −1 C ant contour (Figure S17), where Ω Ar changes are greatest (Figure S3 in Supporting Information S1).Carbonate system nonlinearities work to significantly increase pCO 2 and [H + ] changes below the surface, accounting for >70% of the total changes in waters laden with respiratory byproducts above the 10 μmol kg −1 C ant contour, where changes in these parameters are greatest (Figures S3, S6, and S17 in Supporting Information S1).For the RF, nonlinear contributions to the OA driven changes are predominantly small (<35% of the total RF change), excluding a large dampening signal in the mesopelagic North Pacific that is entirely explained by carbonate system nonlinearities (Figure S17 in Supporting Information S1).In this region the preindustrial DIC:TA ratio is already greater than 1 and increases in DIC from C ant accumulation lead to declines in the RF (Figure S16 in Supporting Information S1).However, adding C ant to the preformed PI values, for which the DIC:TA ratio is much less than 1, causes the RF to increase.Differencing these changes (Equation 4) leads to the large nonlinear reduction in the RF.In summary, carbonate system nonlinearities resulting from natural and anthropogenic carbon pool interactions mitigate C ant -driven changes to Ω Ar and exacerbate changes to pH, pCO 2 , and [H + ] in coherent spatial patterns.The RF responds similarly to pCO 2 in regions where the DIC:TA ratio remains <1 (Figures S3 and S16 in Supporting Information S1).

Potential Implications for Habitat Compression
The largest [H + ] and pCO 2 changes resulting from accumulated C ant in the open ocean have occurred below the sea surface (Figure 2) in mesopelagic waters where organisms are also experiencing reductions in carbonate mineral saturation (Jiang et al., 2015) and dissolved oxygen levels (hypoxia; O 2 ≤ 60 μmol kg −1 ; Breitburg et al., 2018;Oschlies et al., 2018).Elevated in situ pCO 2 has been linked to negative impacts for a suite of sessile and low-vagility invertebrates (Vargas et al., 2022), and there is mixed evidence for adverse impacts on marine fishes through several biological mechanisms (Clements & Hunt, 2015;Esbaugh, 2018;Heuer & Grosell, 2014;Sundin, 2023).One of these mechanisms is a reduction in the efficiency at which CO 2 is expelled during respiration, which can lead to a buildup of CO 2 in the blood (i.e., hypercapnia; Nilsson et al., 2012;Perry & Gilmour, 2006).1,000 μatm is commonly used as a pCO 2 threshold for when respiring marine organisms may experience stress (Arroyo et al., 2022;Feely et al., 2018;McNeil & Sasse, 2016).
The interior ocean pCO 2 we have presented until now reflects the partial pressure that CO 2 would have if it were an ideal gas that was unaffected by the overlying hydrostatic pressure.This is a useful parameter for considering the air-sea pCO 2 gradient that would occur if subsurface waters were instantaneously and quasi-adiabatically brought to the sea surface where pCO 2 and the fugacity of CO 2 (fCO 2 ) are quite similar (Humphreys et al., 2022).However, in the ocean interior, the influence of overlying pressure and non-ideal nature of CO 2 may be an important consideration (Figure S18 in Supporting Information S1).Thus, we compute in situ fCO 2 values (Text S4) and evaluate hypercapnia (hereafter equating to fCO 2 ≥ 1,000 μatm) in the ocean interior and its relation to patterns of low oxygen and Ω Ar .
Hypercapnic conditions can be found at depths shallower than 500 m throughout large expanses of the Pacific Ocean as well as the Bay of Bengal, Arabian Sea, and off the west coast of Africa (Figure 4a).C ant -driven fCO 2 increases in the ocean interior have contributed to this signal, causing widespread shoaling of the hypercapnia horizon (∼180 m along the Pacific transect; green vs. purple line in Figure 4c).A significant volume of hypercapnic water in the Eastern Equatorial Pacific, North Pacific Ocean, and North Indian Ocean is also hypoxic (Figure 4b).In the North Pacific, the hypoxia and hypercapnia horizons closely align with the aragonite saturation horizon at depths ≤200 m (Figure 4d).

Amplified Subsurface Ocean Acidification
OA caused by the air-to-sea transfer of C ant is intuitively thought to be a surface-intensified process, but this is not the case across broad ocean regions and for multiple OA metrics, including pCO 2 , RF, and  values that exceed the upper 50 m mean  ΔCO  ant 2 value (Figure 2) are therefore outpacing the atmospheric pCO 2 perturbation, more than doubling it in some regions of the North Pacific and North Indian Oceans.This is particularly notable since these waters reside below the local maximum annual mixed layer depth and have not been in contact with the atmosphere for some time (Figure S13 in Supporting Information S1).
Aragonite saturation state does not exhibit a subsurface intensified response to C ant accumulation while pH exhibits an inconsistent subsurface intensified response to C ant accumulation in the GLODAPv2.2016bmapped product.There are patches of slightly elevated C ant -driven pH changes below the surface that exceed our calculation uncertainties, and the model analysis (Text S2) indicates that carbonate system nonlinearities are working to induce widespread subsurface pH change amplification (Figure S11 in Supporting Information S1).pH changes reflect relative [H + ] changes (Fassbender et al., 2021;Kwiatkowski & Orr, 2018), so small positive biases in the 2002 subsurface [H + ] values or negative biases in the C ant estimates (which would cause negative biases in  Δ [ H + ]   ) could be muting the large-scale  ΔpH   signal and causing the small and localized subsurface  ΔpH   minima that we find.Alternatively, these patches could reflect unaccounted for uncertainties associated with mapping C ant and/or pH to a global scale to create the GLODAPv2.2016bmapped product.We find that carbonate system nonlinearities are working to amplify the subsurface responses of pCO 2 , RF, [H + ], and pH to C ant accumulation, dominating the overall responses of pCO 2 and [H + ] as well as the response of RF in some ocean regions.Continued OA may cause the patchy  ΔpH  ant to expand into coherent subsurface signals like those found in a recent analysis of 21st century Earth system model projections (Kwiatkowski et al., 2020).Nonlinear carbonate chemistry effects weakly mitigate reductions in Ω Ar , which has a surface intensified response to C ant accumulation.Our findings therefore do not add significant new context to C ant -driven Ω Ar changes relative to previous interior ocean analyses (e.g., Feely et al., 2004;Negrete-García et al., 2019).

Broader Implications
The OA metrics evaluated herein are not routinely incorporated into marine ecosystem models, which are forced with output from Earth system model projections to estimate future ocean biomass changes (Doney et al., 2020;Lotze et al., 2019;Tittensor et al., 2021).However, there is growing interest to understand how climate change, including OA, is impacting the marine environment, and to use this understanding to support ecosystem-based fishery management (Edited by In Bianchi & Skjoldal, 2008;Fletcher et al., 2010;Garcia & Cochrane, 2005;Marshall et al., 2019;Pikitch et al., 2004).A large body of experimental OA work suggests that species from more variable marine environments tend to have more phenotypic plasticity and are thus more tolerant of environmental variability (Boyd et al., 2016).However, most biological OA sensitivity studies use modern or preindustrial near-surface ocean conditions as a baseline mean state from which to conduct their experimental treatments (e.g., Cornwall & Hurd, 2015).Our findings suggest that, in some regions, OA metric changes in subsurface waters with weaker buffering can be significantly larger than surface changes (which may exhibit similar relative changes; Figure S19 in Supporting Information S1), indicating that experimental conditions thought to represent the distant future for the surface ocean may also represent a not-so-distant future for deeper layers of the water column.In addition to mean state considerations, Arroyo et al. (2022) determined that North Pacific waters with amplified subsurface OA changes regularly bathe the continental shelf during upwelling events in the California Current Large Marine Ecosystem (CCLME).Thus, the seasonal range of OA metric extremes on the shelf of the CCLME may have increased by significantly more than in overlying surface waters since the preindustrial period.Inverting the perception that OA impacts attenuate with depth may provide useful new context when designing OA laboratory experiments and interpreting relationships between organismal plasticity and environmental variability (e.g., Vargas et al., 2022).
The nonlinear amplification of interior ocean pCO 2 changes associated with C ant accumulation may have implications for anthropogenic carbon storage in the ocean.C ant enters the ocean due to persistent growth in atmospheric pCO 2 (presently ∼2.5 μatm yr −1 ) that in turn causes persistent growth in surface ocean pCO 2 at a rate that largely tracks the atmosphere; therefore, changes in sea-air pCO 2 disequilibrium (ΔpCO 2 Sea-Air ) have occurred gradually (McKinley et al., 2017).The re-emergence of subsurface waters that have experienced pCO 2 changes >100% larger than global surface ocean changes (Figure S5 in Supporting Information S1) could greatly increase the rate of CO 2 evasion during water mass ventilation.In areas where water masses do not fully equilibrate with the atmosphere prior to subduction (e.g., the Southern Ocean), changes in the rate of CO 2 evasion from the ocean could influence the level of sea-air equilibration reached and thus the amount and efficiency of carbon stored (Eggleston & Galbraith, 2018;Ito & Follows, 2013).With international goals to reduce atmospheric CO 2 levels and eventually achieve net-zero CO 2 emissions (International Energy Agency, 2021), this raises the question: where and when will waters experiencing amplified subsurface pCO 2 changes re-emerge, and could this meaningfully impact ocean carbon storage efficiency?While we have focused on the chemical changes associated with OA, the same concepts hold for other scenarios in which interior ocean carbon or alkalinity gradients are modified; as is the case for many mCDR strategies (National Academies of Sciences, Engineering, and Medicine, 2022).The efficacy and environmental implications of mCDR strategies will be depth-dependent for some key OA metrics.Like the RF, the Alkalinity sensitivity Factor (AF; Takahashi et al., 1993) describes the relative change in pCO 2 associated with a relative change in TA, assuming all other variables are held constant.The AF is elevated at depth so that pCO 2 is more sensitive to alkalinity addition below the surface (Figure S20).This enhanced sensitivity coupled with natural calcium carbonate dissolution (Figure S7) helps to partially compensate for pCO 2 increases resulting from the remineralization of organic matter; however, there is much more carbon added through remineralization than calcium carbonate dissolution (Figure S6).While surface ocean alkalinity enhancement has largely been considered an mCDR strategy to mitigate climate change with the potential co-benefit of OA mitigation, the intentional addition of alkalinity below the surface could mitigate OA impacts that may be more deleterious than previously realized.As these waters will eventually return to the sea surface, mitigating extreme subsurface OA impacts could preemptively mitigate some of the most extreme surface OA impacts as well.Subsurface alkalinity addition could also mitigate future ocean-based CO 2 emissions in regions where subsurface waters with large, C ant -induced pCO 2 perturbations are expected to re-emerge at the sea surface, such as eastern boundary upwelling regions or obduction hot spots.Some of the most abundant mineral sources of TA (e.g., carbonate minerals; Caserini et al., 2022) being considered for use in mCDR do not readily dissolve in modern surface ocean conditions (Kheshgi, 1995) but could dissolve while sinking through undersaturated portions of the water column (Figure 4).

Figure 1 .
Figure 1.Average (a) Revelle sensitivity Factor (RF) and (b) dissolved inorganic carbon (DIC) concentration in the upper 50 m of the ocean.Interior ocean (c) RF and (d) DIC concentrations along 150.5°W (central white line in panels (a) and (b)).Panels a and c as well as b and d use consistent color scales to emphasize the larger interior ocean versus near-surface gradients in these parameters.Data are from the GLODAPv2.2016bmapped product, which is an annual climatology normalized to the year 2002 (Lauvset et al., 2016).White patches in panels (a) and (b) reflect regions lacking data.Panels (a) and (b) also show meridional transects (white lines) in the Atlantic (25.5°W) and Indian (90.5°E)Oceans that will be referenced in subsequent figures.

Figure 2 .FASSBENDER
Figure 2. C ant -induced ocean acidification metric changes (Δ), from the preindustrial period to the year 2002, relative to the local mean change in the upper 50 m for (a-c) pH, (d-f) aragonite saturation state (Ω Ar ), (g-i) partial pressure of carbon dioxide gas (pCO 2 ; μatm), (j-l) [H + ] (nmol kg −1 ), and the (m-o) Revelle sensitivity Factor.Results for the Pacific (150.5°W),Atlantic (25.5°W), and Indian (90.5°E)Oceans with meridional transect locations shown in Figure 1a.Black contours represent C ant (μmol kg −1 ) in the year 2002.Red coloring indicates a subsurface change that is larger in magnitude than the local mean change in the upper 50 m.Stippling indicates where the magnitude of the change is smaller than the uncertainty.

Figure 3 .FASSBENDER
Figure 3. Nonlinear component of ocean acidification (OA) metric changes (Δ) caused by interactions between C ant and natural carbon accumulated from organic matter remineralization and calcium carbonate dissolution.Red and blue coloring indicate where the OA metric responses to C ant are amplified or dampened, respectively, by carbonate system nonlinearities.Results are shown for (a-c) pH (d-f) aragonite saturation state (Ω Ar ), (g-i) partial pressure of carbon dioxide gas (pCO 2 ; μatm), (j-l) [H + ] (nmol kg −1 ), and the (m-o) Revelle sensitivity Factor.Black contours represent C ant (μmol kg −1 ) in the year 2002.Results for the Pacific (150.5°W),Atlantic (25.5°W), and Indian (90.5°E)Oceans with meridional transect locations shown in Figure 1a.Stippling indicates where the magnitude of the change is smaller than the uncertainty.
these parameters exhibited their largest responses to C ant accumulation in the ocean interior, where waters are weakly buffered due to the buildup of remineralized organic material, and C ant concentrations are moderate.Annual mean surface ocean pCO 2 growth largely tracks atmospheric pCO 2 growth, globally and over multi-decadal timescales (McKinley et al., 2017), increasing by ∼90 μatm from the PI to the year 2002 (Figure S1b in Supporting Information S1).Subsurface  ΔCO   2

Figure 4 .
Figure 4. Year 2002 (a) depth of the hypercapnia horizon (fCO 2 = 1,000 μatm) and (b) thickness of overlap in simultaneously hypercapnic and hypoxic (O 2 ≤ 60 μmol kg −1 ) waters.(c) Year 2002 fCO 2 and (d) dissolved oxygen concentration along North Pacific Section 150.5°W.Purple and green lines show the year 2002 and PI hypercapnia horizons (fCO 2 ), respectively.Cyan and yellow lines show the year 2002 and PI aragonite saturation horizons (Ω Ar ), respectively.White lines show the year 2002 hypoxia horizon (O 2 ).

2.1. OA Metric Changes Caused by C ant To
Olsen et al., 2016) of C ant on each OA metric (i.e., pH, Ω Ar , pCO 2 , [H + ], and RF), we rely on a quality-controlled and internally consistent hydrographic data set curated by the Global Data Analysis Project (GLODAP;Olsen et al., 2016)that was previously mapped to a global 1° × 1° horizontal grid with 33 stand- Fassbender et al., 2017)n S1; Equation 1;Fassbender et al., 2017).Similarly large subsurface  ΔCO the North Indian Ocean where C ant penetration is shallow but the vertical pCO 2 :DIC gradient is also shallow and steep (FigureS14in Supporting Information S1).In the North Atlantic and Southern Oceans, where C ant reaches deep into the interior, weaker vertical carbon gradients in the more recently ventilated waters cause subsurface  ΔCO signals in this region to fall within the uncertainties of the calculation.Like pCO 2 and [H + ], the RF exhibits amplified subsurface C ant -induced changes relative to the local mean changes within the upper 50 m, with the largest changes occurring at depths with moderate to low C ant concentrations (Figure2).Poleward of ∼45°N in the Pacific Ocean, there is a large discrepancy in the vertical extent of subsurface   to be negative (Egleston et al., 2010; Fassbender FASSBENDER ET AL. 10.1029/2023GB007843 6 of 14 et al., 2017).The maximum RF is found where DIC approximately equals TA (Figure database (www.glodap.info)by securing funding, dedicating their time to collect and share data, and assembling and quality-controlling those data.AJF was supported by NOAA's Pacific Marine Environmental Laboratory (PMEL).BRC, JDS, and HF were funded through the Cooperative Institute for Climate, Ocean, and Ecosystem Studies (CICOES) under NOAA Cooperative Agreement NA20OAR4320271 and supported by NOAA's PMEL.BRC thanks Kathy Tedesco of NOAA's Global Ocean Monitoring and Observations Program for supporting the Carbon Data Management and Synthesis Grant, which funded his contributions to this project (Fund Ref No. 100007298).MCA was supported by the Eugene Cota-Robles Fellowship administered by UCSC and the NOAA OA Program (Grant NA19OAR0170357).YH was supported by the National Science Foundation (Award No. 2032754) and NOAA PMEL.This is PMEL contribution 5504 and CICOES contribution 2023-1267.