Processes Driving Global Interior Ocean pH Distribution

Ocean acidification evolves on the background of a natural ocean pH gradient that is the result of the interplay between ocean mixing, biological production and remineralization, calcium carbonate cycling, and temperature and pressure changes across the water column. While previous studies have analyzed these processes and their impacts on ocean carbonate chemistry, none have attempted to quantify their impacts on interior ocean pH globally. Here we evaluate how anthropogenic changes and natural processes collectively act on ocean pH, and how these processes set the vulnerability of regions to future changes in ocean acidification. We use the mapped data product from the Global Ocean Data Analysis Project version 2, a novel method to estimate preformed total alkalinity based on a combination of a total matrix intercomparison and locally interpolated regressions, and a comprehensive uncertainty analysis. We find that the largest contribution to the interior ocean pH gradient comes from organic matter remineralization, with CaCO3 cycling being the second most important process. The estimates of the impact of anthropogenic CO2 changes on pH reaffirm the large and well‐understood anthropogenic impact on pH in the surface ocean, and put it in the context of the natural pH gradient in the interior ocean. We also show that in the depth layer 500–1,500 m natural processes enhance ocean acidification by on average 28 ± 15%, but with large regional gradients.


Introduction
Ocean pH is decreasing as a direct consequence of ocean uptake of carbon dioxide (CO 2 ) emitted during the burning of fossil fuels, cement production, and land use changes. This is commonly referred to as ocean acidification (Caldeira & Wickett, 2003;Feely et al., 2004Feely et al., , 2009). The decrease is well documented globally (e.g., Bates et al., 2014;Lauvset et al., 2015) and there is active research focused on understanding its consequences (e.g., Gattuso et al., 2015;Kroeker et al., 2013). The impacts of ocean acidification on marine life have proven to be complex, and there is now ample evidence that other factors such as food availability, predator-prey dynamics, and temperature changes affect an organism's response to ocean acidification (e.g., Kroeker et al., 2013). In a recent metaanalysis, Hurd and Cornwall (2015) highlighted the importance of considering natural gradients in carbonate chemistry for experimental design, in addition to ocean acidification-related chemistry changes.
Understanding the natural processes that determine ocean pH is also important when communicating ocean acidification research. Ocean acidification happens on top of several natural processes, and the combination of the two make some regions particularly vulnerable. One such region is the Southern Ocean, where upwelling of waters rich in remineralized carbon create a local minimum in carbonate ion concentrations in the thermocline. As shown by Negrete-García et al. (2019) this natural gradient means that a very small anthropogenic perturbation is sufficient to shift the aragonite saturation horizon by several hundred meters from one year to the next. A similar effect has been shown by Hauri et al. (2013), Feely et al. (2018), and Franco et al. (2018) for Eastern Boundary Upwelling Systems, which are also regions where natural processes act to enhance the impact of ocean acidification. Such research highlights the fact that ocean acidification decreases natural pH distributions that result from ocean mixing, biological production and remineralization, mineral dissolution, temperature changes, and gas exchange.
Several studies have addressed how natural processes determine the ocean mean distribution of pH (e.g., Feely et al., 2004Feely et al., , 2009), but they are either descriptive or regionally focused. There is a lack of studies quantifying how natural processes control pH globally. Here we assess the subsurface patterns in global climatological ocean pH, and determine the degree to which the pH distribution and its spatial variability are governed by anthropogenic changes, the advection of surface waters into the interior (the so-called preformed component), organic matter remineralization, and carbonate mineral (mostly CaCO 3 , but also, e.g., MgCO 3 ) dissolution. The goal is to properly evaluate how anthropogenic changes and natural gradients collectively act on ocean pH and how this sets the regional vulnerability to future changes in ocean acidification.

Data
Unless otherwise noted we use the Global Ocean Data Analysis Project version 2 (GLODAPv2) mapped climatologies  throughout this study, and all additional calculations are based on this product. The baseline value from which the various pH perturbations are estimated is pH 2002 , which is calculated from the observed total alkalinity (TA) and dissolved inorganic carbon (DIC) in GLODAPv2 , where DIC is first adjusted to year 2002 using the atmospheric perturbation method as detailed in Appendix B of Lauvset et al. (2016), and the assumption that TA does not change over time.
Since few ship-based observations of ocean carbonate chemistry are made in winter there is a recognized seasonal bias in the GLODAPv2 data product , and hence in all calculations made using these data, which makes the data product challenging to use for the surface ocean. No attempt has been made to correct for this bias in this present work, although we note that this bias is primarily an issue for surface measurements and this analysis focuses on the ocean interior. We nevertheless note that previous studies (e.g., Fassbender et al., 2018 ;Landschützer et al., 2018 ;Lenton et al., 2012) have indicated that there is a nonnegligible difference between summertime and wintertime trends in surface ocean pCO 2 , and it is possible that this is also the case for pH. Therefore, the anthropogenic CO 2 impacts on pH (and Ω) estimated here, which align with previous estimates (e.g., Caldeira & Wickett, 2003), should be considered representative of the summertime. Relatedly, warming of the ocean also affects the amplitude of the seasonal cycling and ocean stratification, which has been shown to affect ocean pCO 2 Landschützer et al., 2018). The same effect is likely for ocean acidification, but the GLODAPv2 data product, due to the seasonal bias, cannot be used to assess this effect. This effect has therefore not been accounted for or discussed in detail in this study.
All GLODAPv2 data are available for download through www.glodap.info. In this study, pH and pH changes are always on the total hydrogen scale and, unless otherwise stated, always calculated for in situ temperature and pressure.

Preformed Values
To estimate preformed total alkalinity (TA 0 ), we used a novel combination of the total matrix intercomparison (Gebbie & Huybers, 2011) and the Locally Interpolated Alkalinity Regression (LIARv2) methods (Carter et al., 2018). The approach used here is an evolution of the TA 0 estimate employed by Feely et al. (2002), which relied on basin-specific surface ocean regressions of TA with properties that are conserved in the interior ocean. It addresses two problems with the earlier approach: first, ocean interior ventilation tends to occur during the wintertime when seawater properties are substantially different (e.g., colder) than in the summer. Since the data used to train the regressions have mostly been collected during summer months, this can substantially bias the regressions toward that season, that is, conditions during winter, the time of water mass formation, is not well represented (Carter et al., 2014). Secondly, available TA regressions are valid for specific regions of the surface ocean, and since interior ocean water often originates from distant regions or other ocean basins entirely, selecting which regression to use is problematic. This new approach addresses these challenges by estimating seawater properties just below the shallow summer mixed layers and then propagating these estimates into the global ocean using the total matrix intercomparison approach (Gebbie & Huybers, 2011).
The total matrix intercomparison approach estimates-for a water parcel anywhere in the ocean interiorthe fractions (f i ) of water from each of n total matrix intercomparison grid locations on the ocean surface that comprise the interior ocean mixture. The Locally Interpolated Regressions approach allows estimates that are specific to water masses just below each of the surface locations. Specifically, we use gridded interior ocean potential temperature (θ) and salinity (S) to estimate TA at 100-m depth at each of the n surface locations. Using the 100-m depth surface instead of the 0-m surface limits the impacts of seasonality on the measurements used to train the regressions (Pardo et al., 2011). These estimates now represent the characteristics of water ventilated from each of the n locations, or the TA 0 of each component water mass that collectively comprise the interior mixture. For each ocean interior location on our 1°× 1°× 33 depth levels grid the n individual TA 0 estimates are weighted by f i and then summed, to estimate the average TA 0 of the interior water parcel: To visualize our new preformed alkalinity estimates (TA 0 ) we compare these values to the mapped TA in the GLODAPv2 data product Figures S12 and S13). On the surface the difference is, as expected, very small (1 ± 20 μmol kg −1 , where the uncertainty is 1 standard deviation of the mean; Figure S11a) while it increases with depth to an average 103 ± 6 μmol kg −1 at 3,000 m in the North Pacific Ocean (Figures S12c and S14; uncertainty is 1 standard deviation of the mean) and water mass age ( Figure S1), as a consequence of carbonate mineral dissolution.

Anthropogenic Carbon and Preindustrial Values
Anthropogenic carbon (C anth ) is calculated using the transit time distribution (TTD) method (Hall et al., 2002;Waugh et al., 2006) on all CFC-12 observations in GLODAPv2, following the method used by Lauvset et al. (2016), under the assumption that the ratio between the mean age and width of the TTD is unity (Γ/Δ=1). We used the TA 0 estimated as described in section 2.2 to calculate C anth but note that the TTD based C anth estimates are rather insensitive to how TA 0 is estimated (not shown). C anth was normalized to year 2002 following the method outlined in Lauvset et al. (2016). Preindustrial pH (pH PI ), and aragonite (Ω Ar PI ) and calcite saturation state (Ω Ca PI ) are calculated in CO2SYS (Lewis & Wallace, 1998) version 2.0 for MatLab® (Orr et al., 2018;van Heuven et al., 2009) from DIC PI and TA using the carbonate dissociation constants of Lueker et al. (2000), the bisulfate ion dissociation constant of Dickson (1990), the total boratesalinity relationship of Uppström (1974), and the mapped phosphate and silicate fields from GLODAPv2.
The TTD method used to estimate C anth for this study is known to overestimate C anth in the Southern Ocean due to the assumption of constant equilibrium (He et al., 2018;Khatiwala et al., 2013;Waugh et al., 2006). Waugh et al. (2006) suggests a 20% global overestimate in C anth calculated using the TTD method, while Vazquez-Rodriguez et al. (2009) estimated that the TTD method overestimated C anth in the Southern Ocean by 2-3 μmol kg −1 relative to the Sabine et al. (2004) inventory. As a result, the numbers presented here are likely slightly overestimated.

Decomposition of pH
In the ocean interior pH changes are (i) thermodynamically driven, primarily through changes in temperature (ΔpH T ) and pressure (ΔpH press ), and biogeochemically driven through (ii) net remineralization of organic carbon to inorganic carbon (ΔpH org ) and (iii) net dissolution of CaCO 3 minerals (ΔpH CaCO3 ). The two latter processes do not directly change pH but rather add DIC or release TA (for carbonate dissolution) or titrate it away (for organic matter remineralization). These chemical changes then cause pH to change. In this paper, we consider the measured pH (pH 2002 ) to be the baseline and explicitly calculate how much the interior ocean pH is changed by each process, as well as the accumulation of anthropogenic carbon (ΔpH anth ). Our analysis is based on well-established methods that have been used to decompose DIC and TA in several previous studies (e.g., Cameron et al., 2005 ;Couldrey et al., 2019 ;Eggleston & Galbraith, 2018 ;Gruber et al., 1996 ;Sarmiento & Gruber, 2006;Schmittner et al., 2013). In addition, a more simplified way to decompose interior ocean pH is detailed in Millero (2013), following Park (1969). All components were calculated in CO2SYS (Lewis & Wallace, 1998), using version 2.0 of CO2SYS for Matlab® (Orr et al., 2018;van Heuven et al., 2009), which includes routines for error propagation (Orr et al., 2018) used to estimate uncertainties for each component. The choice of dissociation constants is as described in section 2.3. Our results are visualized as cross sections for the four transects shown in Figure 1: (i) a line through the Nordic Seas along the zero meridian south to the Iceland-Shetland ridge; (ii) the Global Ocean Ship-based Hydrographic Program (GO-SHIP) lines A16 southward in the Atlantic Ocean, SR04 and S04P westward in the Southern Ocean, and P16 northward in the Pacific Ocean; (iii) the GO-SHIP line I09 southward in the Indian Ocean; and (iv) a line that crosses the Arctic Ocean from the Bering Strait to North Pole along 175°W and from the North Pole to the Fram Strait along 5°E. In addition, results are shown as global maps on three different depth horizons (20, 1,000, and 3,000 m).
Preformed preindustrial pH (pH PI,0 ; at in situ temperature and surface pressure) is here interpreted as the preindustrial pH at time of water mass formation. Thus, it is estimated as pH after the impact of the above-mentioned processes have been removed from pH 2002 : where pH 2002 is the pH from GLODAPv2 normalized to year 2002 (see section 2.1), ΔpH CaCO3 represents the change in pH due to carbonate mineral cycling (equation (3)), ΔpH org represents the change in pH due to organic matter remineralization (equation (4)), ΔpH press represents change in pH due to changes in pressure (equation (5)), and ΔpH anth represents the change in pH due to anthropogenic carbon accumulation (equation (6)).
where pH f(DIC,TA) is pH 2002 . This calculation is identical to the used in the GLODAPv2 data product and is repeated in this study only to get the associated uncertainty estimate (see section 2.5). The added subscript P = 0 indicates that the calculation is performed at 0-dbar pressure. Quantifying the impacts of organic matter and carbonate mineral cycling require more steps. First, DIC org and TA org are defined as the dissolved inorganic carbon and total alkalinity, respectively, accumulated at a location through the action of the soft tissue pump, or the remineralization of organic matter: where the stoichiometric ratios are those of Anderson and Sarmiento (1994); the 1.36 coefficient is from Wolf-Gladrow et al. (2007); and AOU, calculated using the Garcia and Gordon (1992) algorithm, is the apparent oxygen utilization, or the difference between the observed dissolved oxygen concentration and the oxygen concentration expected at gas exchange equilibrium (i.e., O 2 sat − O 2 obs ). Second, DIC CaCO3 and TA CaCO3 are defined as the dissolved inorganic carbon and total alkalinity, respectively, accumulated at a location through the action of the hard tissue pump, that is, the net dissolution of carbonate minerals: Third, TA 0 is estimated according to equation (1) while finally DIC 0 is calculated using equation (11): Preformed preindustrial pH (at in situ temperature and surface pressure) can be now be calculated either as pH 2002 adjusted to remove the impacts of the processes considered (equation (2)) or as a function of the preformed preindustrial dissolved inorganic carbon and preformed total alkalinity: Note that ΔpH org is calculated using DIC and TA adjusted for the amounts gained through CaCO 3 cycling (equation (3)) following Sarmiento and Gruber (2006) under the assumption that CaCO 3 dissolution happens after organic matter remineralization has changed both DIC and TA. This minimizes biases brought about by the nonlinear dependence of pH to DIC and/or TA perturbations.
As mentioned, all calculations are conducted at in situ temperature, under the assumption that this is conserved as water circulates through the interior ocean. Effects of hydrothermal heating and adiabatic expansion and compression of seawater are small and neglected. For example, the impact on pH of the heating water masses experience as they are transported from surface to higher pressures (adiabatic heating), is very small. Typically, a change of pressure from surface to 4,000 m generates a heating of about 0.4°C, which would lead to a pH drop of~0.006 units (as the pH change per 1°C is~0.015 under constant DIC and TA). The impacts of ocean mixing are captured in our preformed pH distributions (at in situ temperature). However, we caution that temperature plays a major role in driving the processes that lead to surface ocean variability in carbonate chemistry (Jiang et al., 2015) while only having a small net impact on the pH of surface waters that are well-equilibrated with the atmosphere (Jiang et al., 2019).

Uncertainty Analysis
The uncertainties in each input variable and parameter are listed in Table 1, and the error propagation equations for equations (3)-13 are given in equations S2-S15 in the supporting information. Following recommendations by Orr et al. (2018) these uncertainties are defined as standard uncertainties (1σ), and we assume that the uncertainties in DIC and TA are uncorrelated. Since we include the mapping error the assumption of uncorrelated uncertainties is not strictly true since where mapping error for DIC is large it is also large for TA. This leads to an overestimation of the overall uncertainty (Orr et al., 2018).
Throughout the text and in all figures the final combined uncertainties ( Figures S4-S8) are presented as 2σ. Some important caveats regarding our uncertainty analysis are detailed in the supporting information.

pH Gradients Due to Natural Processes
Interior ocean pH 2002 gradients are substantially larger than surface pH 2002 gradients, with an in situ pH 2002 range between 7.52 ± 0.05 and 8.33 ± 0.10 in the interior ocean (Figures 2a and S4). Mostly, this range is due to biological processes that affect the interior ocean pH as outlined in section 2.4, and pressure effects. Pressure has a significant impact on pH in the deep ocean interior, primarily due to the pressure effects on the dissociation constants for carbonic acid (K 1 ) and bicarbonate ion (K 2 ), and causes a decrease in pH with depth ( Figure 3a). Along our section (Figure 1), which is never deeper than 5,500 m, ΔpH press has a maximum impact of −0.27 ± 0.11, that is, the chemical speciation changes that result from increased pressure at depth decrease the pH in the deepest parts of our section by 0.27.
Organic matter remineralization leads to a decrease in interior ocean pH (Figure 3b) because remineralization increases the concentration of inorganic carbon. Of the four components that create the interior ocean pH gradients (equation (2)) ΔpH org has the largest contribution overall, and has its largest impacts in the northern Pacific and the northern Indian Oceans where subthermocline water masses are typically older than 1,000 years ( Figure S1): ΔpH org is up to −0.79 ± 0.28 and −0.61 ± 0.23 in these regions, respectively ( Figure 3b). As we use AOU to estimate organic matter remineralization there is naturally also a close correlation between maximum ΔpH org and oxygen minimum zones, which are often not caused by high rates of organic matter remineralization, but rather lack of mixing and ventilation due to slow or stagnant circulation (Keeling et al., 2010). Along our transect, the spatial pattern ΔpH org is also very similar to that in Mapping uncertainties are taken from the GLODAPv2 mapped data product . d Estimated based on the results of He et al. (2018). e Assumed based on range of published, commonly used stoichiometric ratios (e.g., Redfield et al., 1963;Takahashi et al., 1985;Li & Peng, 2002;Körtzinger et al., 2001). f 0.3% following Garcia and Gordon (1992). g Estimated based on published values of this coefficient (Kanamori & Ikegami, 1982;Wolf-Gladrow et al., 2007).
h Default uncertainty as given in the errors.m subroutine of the CO2SYS program (Orr et al., 2018

Global Biogeochemical Cycles
LAUVSET ET AL.
both water mass age ( Figure S1) and DIC org ( Figure S9) since the amount of remineralized carbon in seawater increases as the product of the age of the water mass and the average net rate of remineralization. That the age of the seawater is an important factor (r = 0.65) can clearly be seen in the steady decrease in the magnitude of ΔpH org moving from south to north across the Pacific (Figure 3b). Both from the AOU (not shown) and the water mass age, the Southern Ocean (south of 45°S) appears quite recently ventilated, as expected due to active water mass formation in this region. Antarctic Intermediate Water and Subantarctic Mode Water are formed in the Antarctic Circumpolar Current, and spread northward to ventilate the Southern Hemisphere at~1,000-m depth. These water masses show up with comparatively small ΔpH org signals in all ocean basins (Figure 4b). There is nevertheless a significant subsurface ΔpH org component (up to −0.37 ± 0.15) due to upwelling of quite old Pacific Deep Water and Upper Circumpolar Deep Water (Talley, 2013) near the polar front of the Antarctic Circumpolar Current (Figure 3b).
The importance of the rate of remineralization, approximated by dividing ΔpH org by water mass age (Figures S10 and S11), can also be seen, for example, in the equatorial Atlantic where a ΔpH org maximum is found in northward penetrating Antarctic Intermediate Water, which is younger than the North

Global Biogeochemical Cycles
Atlantic Deep Water below it (Figures 4b and 4c), but shallower and subjected to higher rates of remineralization (Figures S11b and S11c). Local ΔpH org maxima found in the upper 1,000 m of the ocean (Figures 3b and 4a-4c) are all associated with the thermocline or upwelling of older water and highlight that the largest vertical gradients in ΔpH org are in the upper ocean, which is consistent with the Martin curve (Martin et al., 1987). Upwelling deep waters show up as patches of strong remineralization signals at 100 m on both sides of the equator in the Pacific, and along the U.S. western seaboard (Figure 4a). The deep ocean, on the other hand, generally has a very small, and spatially uniform, rates of remineralization ( Figure S11c), but also here there are regional differences with well-ventilated areas generally having larger rates than old water masses. Related to this there is a clear inverse correlation (r = −0.69) between the rate of remineralization and C anth (not shown).
In contrast to organic matter remineralization, CaCO 3 dissolution leads to an increase in pH (Figure 3c) because dissolution of CaCO 3 releases carbonate ions, thereby elevating TA more than DIC and increasing the buffer capacity (Feely et al., 2002(Feely et al., , 2012Zeebe & Wolf-Gladrow, 2001). Note that the spatial pattern in ΔpH CaCO3 closely resembles that of ΔDIC CaCO3 (see Sarmiento & Gruber, 2006, Figure 9.3.4), although  Figure 1, of (a) the change in pH incurred by pressure changes (ΔpH press ; equation (5)), (b) the change in pH incurred by organic matter remineralization (ΔpH org ; equation (4)), (c) the change in pH incurred by the dissolution of CaCO 3 in the interior ocean (ΔpH CaCO3 ; equation (3) our maximum is somewhat shallower than that shown by Sarmiento and Gruber (2006). The impacts of carbonate mineral cycling are typically much smaller than the impacts of organic matter cycling both because more carbon is cycled as organic matter than as carbonate minerals and because the effects on pH of increased TA (2 per CaCO 3 dissolved) are partially offset by the effects of increased DIC (1 per CaCO 3 dissolved). Interior ocean net carbonate dissolution lead to an increase in pH of at most 0.24 ± 0.12 in the North Pacific, and 0.22 ± 0.10 in the northern Indian Ocean. In both these regions, the maximum impact is centered around 2,500 m. In the North Pacific, the calcite saturation horizon is very shallow (200-700 m), although still approximately 100-300 m deeper than the aragonite saturation horizon (Figure 3c), which is consistent with previous studies (Feely et al., 2002(Feely et al., , 2012. Though wellknown and much described, the shallow and similar depths for aragonite and calcite saturation horizons are unusual since everywhere else the calcite saturation horizon is deeper than 3,000 m, and~2,000 m deeper than the aragonite saturation horizon. This pattern is ultimately governed by the remineralization of organic matter, the DIC accumulation in the oldest water masses leads to greater decreases in pH-as discussed above-and corrosive conditions for calcite and aragonite. The global pattern of ΔpH CaCO3 therefore resembles the pattern of ΔpH org , but with an opposite sign and reduced magnitude for reasons mentioned above. Interestingly, the area of high ΔpH CaCO3 found in the northern Indian Ocean does not have a corresponding shallow calcite saturation horizon, as found in the North Pacific, but the aragonite saturation horizon does shoal toward the surface. Distributions of pH PI,0 P = 0 , the preindustrial pH at time of water mass formation, which are subsequently only changed by mixing in the interior ocean, are displayed in Figures 3d and 4g-4i. pH PI,0 P = 0 is highly homogeneous in the ocean interior with an average of 8.20 ± 0.14 at both 1,000 and 3,000 m (Figures 4h  and 4i). Given the sensitivity of pH to temperature, and the relatively large temperature gradients near the surface (i.e., at 100 m), the surface gradients in pH PI,0 P = 0 are large (Figure 4g). However, below 1,000 m pH PI,0 P = 0 is 8.20 ± 0.13 throughout the transects chosen here (Figure 3d). This indicates that the major natural processes, which add or remove DIC and TA in the ocean interior and thus change pH, have been accounted for and quantified fairly accurately. As an additional control for the validity of our decomposition we use the estimated TA 0 (equation (1)) and DIC PI,0 (equation (11)) to calculate preformed preindustrial partial pressure of CO 2 (pCO 2 PI,0 ; Figure S14). Given that DIC PI,0 accumulates all sources of error, it lends confidence in our preformed properties, and thus our decomposition, that the difference between our pCO 2 PI,0 estimate and the approximate preindustrial atmospheric CO 2 (xCO 2 PI,atm~2 80 ppm) is consistent with theoretical expectations in having higher undersaturation signals in the high latitudes compared to the lower latitudes. Note that nowhere is our pCO 2 PI,0 significantly different from xCO 2 PI,atm .

Anthropogenic Changes
The anthropogenic change in ocean pH is due to the ocean's increasing absorption of CO 2 from the atmosphere. This uptake is driven by the increasing partial pressure of CO 2 in the atmosphere and has led to more carbon being stored in the ocean. The estimates of the impact of anthropogenic CO 2 changes on pH (Figure 5b) reaffirm the large and well-understood anthropogenic impact on pH in the surface ocean (Jiang et al., 2019), and put it in the context of natural pH gradients in the interior ocean. In the top 200 m ΔpH anth is >80% of ΔpH org , and in the convergence zones ΔpH anth is at least 40% of ΔpH org down to 800 m. Along the section shown here, the 0.05 ΔpH anth contour only extend deeper than 1,000 m in the North Atlantic, Arctic, and parts of the Southern Ocean (Figures 2a and 5b)-regions where the water column inventory of C anth (Figure 5a) is large. Thus, 1,000 m is generally the depth below which natural gradients dominate and anthropogenic changes are small.
This paper focuses mainly on interior ocean carbonate chemistry gradients, but it is important to consider surface ocean changes for context and because anthropogenic changes in pH originate at the surface. Surface ocean pH is clearly lower in year 2002 than it was in preindustrial times (Figures 6a and 6b). The decrease in surface waters is a globally quite uniform 0.10 ± 0.10 (Figure 6c and see Text S3 in the supporting information for notes on the uncertainty), although the uncertainty is spatially more variable than the change itself ( Figure S7c). Since the pH scale is logarithmic this change represents a 25% increase in surface ocean acidity (total hydrogen ion concentration) relative to preindustrial times (approximately year 1800

Anthropogenic Change and Its Drivers
Anthropogenic ocean pH changes are a direct function of the accumulation of C anth in the ocean so any quantification of ΔpH anth depends to some degree on how C anth is estimated. In this study we use the  Figure 1, of ΔpH anth (equation (6)). The vertical axis shows depth (m) below sea level.

Global Biogeochemical Cycles
TTD method (see section 2.3) and the estimated water column inventory of C anth (Figure 5a) of 167 ± 29 PgC is (just) within the uncertainties of the estimate by Khatiwala et al. (2013) for 2010 (155 ± 31 PgC) and the more recent of estimate of 160 ± 20 PgC for 2010 by Gruber et al. (2019) based on the eMLR(C*) approach (Clement & Gruber, 2018). Our estimate is for an earlier year (2002) and yet is the highest of the three, but regardless, the comparability with previous estimates, along with the spatial pattern being consistent with previous work (e.g., Khatiwala et al., 2013 ;Sabine et al., 2004), lends confidence to the estimates of the impact of anthropogenic CO 2 on pH presented here.
Ocean pH is changing as a direct response to accumulation of C anth , but the magnitude of pH change for a given change in DIC (∂pH/∂DIC) is dependent on how much DIC and TA are present (Figure 7). Near the surface, ∂pH/∂DIC is smaller than at depth, meaning that a surface change in C anth will lead to a smaller change in pH than a deep change. Figure 8a, comparing the ΔpH anth with the accumulated C anth , shows that the change in pH for any given C anth concentration is larger (−0.0026 ± 0.000013 pH units per μmol kg −1 C anth (where the uncertainty is the 95% confidence interval around the slope)) when C anth is low (<40 μmol kg −1 ), than where C anth is high. This reflects that the DIC/TA ratios, and thus ∂pH/∂DIC, are largest in the deep ocean where the anthropogenic carbon concentrations are low. There is a similar difference when looking only at the surface waters (i.e., where C anth is high) with larger ∂pH/∂DIC in cold water than in warm water (compare Arctic and Tropics in Figure 7). For Ω Ar , on the other hand, there is no apparent breakpoint in the rate of change with C anth (Figure 8b) which is because the sensitivity of a change in Ω Ar for a given change in DIC (∂Ω Ar /∂ DIC) is strongest at the surface and in warm water ( Figure S17). The relationship visualized in Figure 7 highlights that anthropogenic changes happens on top of a natural gradient, and that the size of this gradient will influence ΔpH anth , both on the surface and in the interior ocean. To assess how the natural processes affect the anthropogenic impacts we calculate pH anth no¯nat from the preformed values (equation (14)) and compare that with the outcome of equation (6). The subscript "no nat" here indicates that the natural processes have been removed.

Figure 7.
Visualization of the sensitivity of change in pH for a given change in DIC (∂pH/∂DIC) as a function of DIC and TA. ∂pH/∂DIC is calculated from the GLODAPv2 data using the derivnum function. This function is part of the CO2SYS software, and has been modified here to output pH (rather than H+). Note that the sensitivity has been calculated for surface pressure (0 dbar) and 25°C in order to visualize the effects of DIC and TA only.
The difference between ΔpH anth and pH anth no¯nat is presented in Figure 9. Clearly, the inclusion of DIC from the remineralization of organic matter (DIC org ) in the estimate gives a greater pH change for any given increase in C anth . This result reflects the changing ∂pH/∂DIC as water circulates in the deep ocean and accumulate DIC org : C anth is added to the ocean at the surface, water then sinks out of the surface, DIC org is added as the waters circulate in the deep ocean, ∂pH/∂DIC increases because of the increased DIC, and the increased sensitivity leads to a strengthening of the the initial pH decline brought about by the uptake of C anth at the surface. Note that comparing equations (6) and (14) does not directly tell us how much more sensitive pH is to addition of DIC org , which would have to be quantified by modifying equation (4), only how much the initial acidification is increased when remineralization subsequently happens. Increased TA from carbonate mineral dissolution, the other natural process discussed in this paper, would be expected to decrease ∂pH/∂DIC and thereby reduce the pH change due to addition of C anth , but these chemical changes are more than compensated by the release of CO 2 through organic matter remineralization, which is clearly seen when comparing ∂pH/∂DIC at 3,000 m in an ocean without biology (labelled 3,000 m; TA 0 , DIC 0 in Figure 7) and ∂pH/∂DIC at 3,000 m in an ocean with biology (labelled 3,000 m in Figure 7). Hereafter we refer to the strengthened ΔpH anth due to the addition of DIC org as an enhancement.
We find that in the depth layer 500-1,500 m the net enhancement is on average 28 ± 15% along our section, but there are large regional and vertical gradients (Figure 9). In the deep ocean, except in the Southern Ocean, the differences between pH anth no¯nat and ΔpH anth are generally very small (<0.005; Figure 9a) and have therefore been removed from further analysis (Figure 9b). In the Southern Ocean there is a significant amount of C anth at the bottom, and here the ocean acidification enhancement decreases from 71 ± 4% in the 800-1,500-m range to 25 ± 4% deeper than 4,000 m. In the upper ocean (<1,500 m), however, the largest enhancements are found in low-oxygen regions like the equatorial Atlantic where the difference is 48 ± 13%. In the Indian Ocean the strongest enhancement (70 ± 10%) is found between 1,200 and 1,750 m just north of   (6)) and the same when the natural processes are removed first (equation (14)). (b) Vertical cross section showing the relative difference between equations (6) and (14), that is, (equations (6)- (14))/equation (6), presented as a percent change. Note that everywhere the numerator is less than 0.005 has been removed. This is because the denominator is also generally very small in these areas and this creates misleadingly large enhancement. The vertical axis shows depth (m) below sea level.
the Broken Plateau, where high AOU corresponds with the shoaling of the aragonite saturation horizon. In the Pacific Ocean the maximum along our section is 35 ± 14% at 300-1,750 m. While these relative changes appear large, the largest relative enhancements are found in areas of low oxygen that tend to have low concentrations of C anth . However, as C anth continues to be added to the ocean and propagate into poorly ventilated portions of the ocean this enhancement will grow more pronounced. In this regard the regional patterns are important to quantify as they render some ecosystems more vulnerable than others.

Ecosystem Vulnerabilities
Research is ongoing as to at what point the decrease in pH, Ω Ar , and Ω Ca will exceed ecosystems' natural thresholds for sustainability. In that respect, it is interesting that a recent study (Sulpis et al., 2018) found that there is already an observable anthropogenically induced CaCO 3 dissolution at the ocean bottom in the North Atlantic, where deep-sea corals persist in comparatively high-pH deep waters with little naturally accumulated carbon from organic matter remineralization. By contrast, in the North Pacific Ocean pH is naturally low because of the strong effect of organic matter remineralization which is only partially compensated by CaCO 3 dissolution. The northern North Pacific also has a very shallow aragonite saturation horizon (<200 m in some places). Thus, despite the anthropogenic impact overall being smaller in the North Pacific compared to the North Atlantic, it is possible that Pacific ecosystems are more vulnerable to ocean acidification because the changes compound on naturally low pH; that is, they are already close to CaCO 3 undersaturation. However, assessing the impact of ocean acidification on marine organisms is difficult since the impacts are often species specific and regionally varying (e.g., Fabry et al., 2008 ;Gattuso et al., 2015 ;Kroeker et al., 2013), but what is clear is that heavily calcified organisms are particularly sensitive, while organisms that are able to move themselves have higher tolerances (Kroeker et al., 2013). It has been shown that it is the CaCO 3 saturation state rather than [CO 3 2− ] or pH which controls calcification (e.g., Bednaršek et al., 2017;Fabry et al., 2008;Osborne et al., 2019;Waldbusser et al., 2015), and recent research from the North Pacific indicate that some calcifying species are quite vulnerable to anthropogenic change (Bednaršek et al., 2014(Bednaršek et al., , 2017Feely et al., 2016;Osborne et al., 2019). Other studies (e.g., Fabry et al., 2008) have found that organisms adapted to strong natural gradients in carbonate chemistry have better mechanisms for buffering such variations, but it is still unclear whether this means they are also more resilient to anthropogenic changes that shift the natural gradient toward lower pH and shallower saturation horizons. Furthermore, naturally low-pH, high-DIC waters have a higher ∂pH/∂DIC and are thus more sensitive to further C anth accumulation (Figures 7-9). In the Southern Ocean, a region sensitive to anthropogenic changes due to the natural shallow minimum in (CO 3 2− ), a recent model study (Negrete-García et al., 2019) shows that there is a high likelihood of a very rapid emergence of shallow (<100 m) saturation horizon for aragonite, largely due to the accumulation of CO 2 in the thermocline, which decreases the natural minimum in [CO 3 2− ] enough to cause undersaturation.
While this paper focuses on pH in the interior ocean the anthropogenic change is largest at the surface; thus, any ecosystem impact is likely to first be seen there. As an example application of the quantified anthropogenic changes in pH and Ω presented in this study we have estimated the potential loss of suitable habitat for shallow calcifying corals. In today's ocean coral reef ecosystems do not thrive where the water has Ω Ar < 3 (Guinotte et al., 2003), and Figure 10 shows the extent of the ocean areas with surface Ω Ar ≥ 3 in 2002 ( Figure 10a) and in preindustrial times (Figure 10b). The potential habitat loss, defined as the preindustrial area minus the area in 2002, is 35,300,000 km 2 . The calculation uncertainty for Ω Ar is substantial, largely due to the large uncertainty in the solubility product (Orr et al., 2018). Thus, the smallest possible extent in preindustrial times approximately equals the largest possible extent in 2002 and the uncertainty in the loss of area (ranging from 35,000,000 to 83,000,000 km 2 ) is very large. Regardless of the change in habitat with Ω Ar ≥ 3 an interesting feature in Figure 10 is that there are large regions today which appear to be very vulnerable to further anthropogenic changes. If the minimum 2002 area is used then the South China Sea and the East Indian Archipelago as well as most of the west coast of the Americas all have Ω Ar < 3. If the maximum area is used, these regions mostly have Ω Ar > 3, but the fact that just the uncertainty in the global mapped Ω Ar yield such large differences indicate that these regions are just within the "healthy" zone and thus particularly vulnerable to additional anthropogenic changes. The above estimate of area loss, calculated using the best available observations of ocean carbonate chemistry, highlights the need for reducing the uncertainty in both observations and calculations. Until then, fully quantifying the environmental stress all ocean ecosystems are subject to, both today and in the future, remain difficult.

Caveats
There are some notable caveats regarding possible sources of bias in our study. First, it is clear from much previous work that pH PI,0 P = 0 is sensitive to how DIC org and TA org are calculated. AOU, which is used in our study, is the most commonly used method because it does not require knowledge about preformed values. However, previous studies (e.g., Carter et al., 2014;Ito et al., 2004;Talley et al., 2003;Russell & Dickson, 2003) have raised questions about the appropriateness of AOU as a proxy for organic matter remineralization given the assumption that a water mass is fully saturated with atmospheric oxygen at the time of subduction. This has for a long time been known to be incorrect (e.g., Redfield et al., 1963) since the degree to which saturation is reached depends on many factors including, but not limited to, wind speed, cooling, ice cover, and biological processes on the surface (e.g., Duteil et al., 2012). In addition, the use of AOU as a proxy for organic matter remineralization does not account for the fraction that occurs through denitrification. Deutsch et al. (2001) estimated the impacts of water column denitrification in the North Pacific, finding a maximum nitrate loss in oxygen-deficient zones of~10 μmol kg −1 , which implies that our AOU-based DIC org is off by a comparable number. Denitrification can therefore be assumed to result in a ΔpH org estimate error of~0.03 within oxygen-deficient zones and considerably less elsewhere. Second, constant stoichiometric ratios have been used throughout this study, while there is ample evidence that these ratios vary both spatially and temporally (e.g., Arrigo et al., 1999;Frigstad et al., 2011;Geider & La Roche, 2002;Li & Peng, 2002). However, as discussed in Deutsch and Weber (2012), synthesis efforts show that there appears to be a global mean stoichiometry similar to the classical Redfield ratio. Using spatially varying stoichiometry would likely change the results presented here, so we conservatively assume a 30% uncertainty

10.1029/2019GB006229
Global Biogeochemical Cycles applied to our globally uniform r C:O value (Table 1) to account for the possible latitudinal range of variability in stoichiometric ratios.

Summary and Conclusions
We decompose pH into components describing the anthropogenic impact (ΔpH anth ), organic matter remineralization (ΔpH org ), CaCO 3 cycling (ΔpH CaCO3 ), thermodynamics (ΔpH T and ΔpH press ), and ocean circulation (pH PI,0 P = 0 ). The ΔpH org component is the most important (maximum reduction of 0.79 ± 0.26) in explaining the interior ocean pH gradients, while ΔpH CaCO3 has a significant effect only at quite large depths (centered around 2,500 m) in the North Pacific and the northern Indian Oceans. The Atlantic Ocean sector of the Southern Ocean, which has a naturally shallow aragonite saturation horizon, is one area where the gradient between the surface ocean pH (~8.2) and the naturally low subsurface pH (~7.9) is quickly vanishing due to invasion of anthropogenic carbon. Here pH is naturally low due to strong organic matter remineralization, and anthropogenic changes are penetrating very deep (ΔpH anth~0 .07 at 1,000 m), and this is an area recently identified as being particularly sensitive to small changes in ocean carbonate chemistry (Negrete-García et al., 2019). When removing the organic matter remineralization and CaCO 3 cycling from the estimate ΔpH anth becomes smaller, indicating an enhancement of anthropogenic pH changes by natural processes, which is likely to continue and grow stronger as ocean acidification continues. With the exception of the Southern Ocean, we currently see this enhancement only in the upper ocean (<1,500 m) and most pronounced in low-oxygen areas. Anthropogenic emissions of CO 2 have led to both decreased pH and decreases in the saturation states of CaCO 3 minerals. However, given the uncertainties in calculating the ocean carbonate chemistry it is difficult to conclude definitively how large this impact is on ocean ecosystems. Thus, fully quantifying the environmental stress ocean ecosystems are subject to both today and in the future through the reduction of uncertainty in the observations and calculations remains a priority.