Interactive comment on “ Natural variability in the surface ocean carbonate ion concentration ”

Abstract. We investigate variability in the surface ocean carbonate ion concentration ([CO32−]) on the basis of a~long control simulation with an Earth System Model. The simulation is run with a prescribed, pre-industrial atmospheric CO2 concentration for 1000 years, permitting investigation of natural [CO32−] variability on interannual to multi-decadal timescales. We find high interannual variability in surface [CO32−] in the tropical Pacific and at the boundaries between the subtropical and subpolar gyres in the Northern Hemisphere, and relatively low interannual variability in the centers of the subtropical gyres and in the Southern Ocean. Statistical analysis of modeled [CO32−] variance and autocorrelation suggests that significant anthropogenic trends in the saturation state of aragonite (Ωaragonite) are already or nearly detectable at the sustained, open-ocean time series sites, whereas several decades of observations are required to detect anthropogenic trends in Ωaragonite in the tropical Pacific, North Pacific, and North Atlantic. The detection timescale for anthropogenic trends in pH is shorter than that for Ωaragonite, due to smaller noise-to-signal ratios and lower autocorrelation in pH. In the tropical Pacific, the leading mode of surface [CO32−] variability is primarily driven by variations in the vertical advection of dissolved inorganic carbon (DIC) in association with El Nino–Southern Oscillation. In the North Pacific, surface [CO32−] variability is caused by circulation-driven variations in surface DIC and strongly correlated with the Pacific Decadal Oscillation, with peak spectral power at 20–30-year periods. North Atlantic [CO32−] variability is also driven by variations in surface DIC, and exhibits weak correlations with both the North Atlantic Oscillation and the Atlantic Multidecadal Oscillation. As the scientific community seeks to detect the anthropogenic influence on ocean carbonate chemistry, these results will aid the interpretation of trends calculated from spatially and temporally sparse observations.


Introduction
The global ocean has absorbed ∼ 30 % of the carbon dioxide (CO 2 ) released by human activities since 1765 (Ciais and Sabine, 2013).While ocean uptake of CO 2 plays a key role in mitigating anthropogenic climate change, it also modifies ocean carbonate chemistry (Feely et al., 2004).The dissolution of excess CO 2 in the surface ocean drives an increase in the dissolved inorganic carbon (DIC) concentration without changing the alkalinity (Alk).The result is a surface ocean characterized by decreasing carbonate ion concentration ([CO 2− 3 ]) and pH (Feely et al., 2009).This acidification of the surface ocean reduces the saturation state of the calcium carbonate minerals calcite and aragonite ( calcite and aragonite , respectively) and may reduce biogenic calcification and enhance calcium carbonate dissolution (Doney et al., 2009).
Observations collected at sustained open ocean time series stations (e.g., HOT, BATS) indicate significant anthropogenic changes in surface DIC, [CO 2− 3 ], pH, and aragonite relative to background natural variability (Bates et al., 2014).The detection of statistically robust trends in carbonate chemistry at these stations benefits from frequent sampling (3-16 times per year), long records (15 to 30 years), and low natural variability (Le Quéré et al., 2000;Brix et al., 2004;Bates et al., 2014).In the rest of the global ocean, however, sparse spatial and temporal sampling, coupled with potentially large natural variability challenges the detection of anthropogenic changes in carbonate chemistry from observa-Published by Copernicus Publications on behalf of the European Geosciences Union.Sutton et al. (2014) report decreasing pH from 1997 to 2011 using mooring observations, but they attribute approximately 40 % of this decrease to natural variability.Based on measurements from repeat hydrographic surveys, Feely et al. (2012) report an average decrease in aragonite and calcite of 0.34 % yr −1 in the Pacific Ocean.In the South Pacific, the trend is primarily driven by uptake of anthropogenic CO 2 , while the trends in the North Pacific Subtropical Gyre and the California Current are attributed to natural variability in ocean circulation.On the global scale, Lauvset et al. (2015) find a mean rate of decrease in surface ocean pH of 0.0018 yr −1 over 1991-2011, using observations of f CO 2 aggregated into 17 biogeographical biomes.They find a substantial amount of interannual variability in pH in many of the biomes (RMSE ranging from 0.01 and 0.04 pH units) that is of a similar magnitude to the cumulative trend in pH.
Internal climate variability arises from the coupled interaction of atmospheric, oceanic, terrestrial, and cryospheric processes (Deser et al., 2012a) and complicates our ability to detect anthropogenically forced trends from sparse observations.In the tropics, the dominant mode of internal climate variability is the El Niño-Southern Oscillation (ENSO).In the extratropics, three major climate modes drive variability: the North Atlantic Oscillation (NAO), the Pacific Decadal Oscillation (PDO), and the Southern Annular Mode (SAM).Studies conducted using ocean physical and biogeochemical models run in hindcast mode (i.e., forced with the historically observed atmospheric state) suggest that ENSO, NAO, PDO, and SAM impact regional ocean biogeochemistry.Le Quéré et al. (2000) and Long et al. (2013) find reduced CO 2 outgassing in the tropical Pacific during El Niño events as a result of changes in dynamics reducing the vertical advection and diffusion of DIC into the surface ocean.In the North Atlantic, the NAO drives shifts of the subpolar/subtropical intergyre boundary that affect the vertical and lateral advection of DIC and air-sea CO 2 flux (Thomas et al., 2008).On the basis of seven biogeochemical models run in hindcast mode, McKinley et al. (2006) show that the positive phase of the PDO is associated with an increased surface DIC tendency in the subtropical and western subpolar gyres of the North Pacific.In the Southern Ocean, multiple hindcast modeling studies find large interannual variability in surface DIC and Alk driven by the SAM (Lenton and Matear, 2007;Lovenduski et al., 2007;Verdy et al., 2007).
Model hindcast studies are useful for quantifying the impact of climate variability on ocean carbonate chemistry, but are limited in their temporal scope to the period of time for which we have abundant observations of the global atmosphere (typically 1948 to the present day).Large-scale modes of climate variability such as PDO, NAO, and SAM have spectral power at low frequencies.While hindcast studies can capture the observed chronology of these modes, they cannot capture the full spectrum of internal variability in the climate system.Long model simulations (order 1000 years) can cap-ture multiple realizations of climate variability on decadal and multi-decadal timescales, and have shown to be useful in the study of ocean carbon cycle variability on these timescales (Doney et al., 2006;Séférian et al., 2013Séférian et al., , 2014;;Keller et al., 2014;Lehner et al., 2015;Resplandy et al., 2015).
Here, we assess the influence of internal climate variability on surface ocean carbonate chemistry by analyzing output from a 1000-year control simulation of a coupled Earth System Model.Interaction between the model's atmosphere, ocean, terrestrial biosphere, and cryosphere generates internal climate variability on timescales ranging from interannual to multi-decadal and longer.We aim to quantify and mechanistically understand the drivers of variability in surface ocean carbonate chemistry on these timescales.In doing so, we will gain perspective on the statistical confidence in the anthropogenic carbonate chemistry trends reported in the literature.
Our study builds upon two recent studies of natural variability in ocean carbonate chemistry from long integrations of Earth System Models (ESMs).Friedrich et al. (2012) use MPI-ESM to quantify the natural variability in surface aragonite and compare it to the anthropogenic trend over the years 800-2099.They suggest that recent anthropogenic trends in surface aragonite exceed natural variability by 30 times on regional scales, but do not focus on detectability in the observational record or on the mechanisms driving variability.Séférian et al. (2013) analyze output from a fully coupled 1000-year control simulation of IPSL-CM5A-LR and describe decadal to multi-decadal variability in air-sea CO 2 flux and its driving factors in the North Atlantic, North Pacific, and the Southern Ocean.They find that a large fraction of the variance in CO 2 flux is driven by internal climate variability in the various regions, due to circulation-mediated variability in the upwelling of DIC to the surface ocean, but only very briefly discuss the implications of this for carbonate ion variability.Here, we analyze output from a 1000-year control simulation of the Community Earth System Model (CESM) with a focus on quantifying and understanding the drivers of variability in surface ocean carbonate chemistry.Unlike the simulation analyzed in Friedrich et al. (2012), the CESM control simulation does not include any external forcing, such as anthropogenic CO 2 emissions, volcanic eruptions, or solar variability.This allows us to cleanly ascribe surface carbonate variability to internal variability (unforced natural variability) of the physical climate system.Further, unlike the simulations analyzed in Friedrich et al. (2012) and Séférian et al. (2013), we analyze output from a simulation with constant, prescribed atmospheric CO 2 concentration, so variability in ocean biogeochemistry is only affected by the physical state of the atmosphere and ocean and not by variability in atmospheric CO 2 .This simplifies our quantification of the mechanisms driving variability in carbonate chemistry.Finally, we focus on variability in surface ocean [CO 2− 3 ], as it is the primary source of variability in aragonite , and is likely to be influenced by internal climate variability ([CO 2− 3 ] ≈ Alk − DIC).As such, our results will be useful for determining the detectability of anthropogenic trends in aragonite over background internal climate variability on a global scale.

Model description
We analyze output from a 1000-year pre-industrial control simulation of the Community Earth System Model.CESM is a state-of-the-art coupled climate model consisting of atmosphere, ocean, land, and sea ice component models (Hurrell et al., 2013).The atmosphere model is the Community Atmosphere Model, version 4 (CAM4), with a horizontal resolution of 1.25 • × 0.9 • and 26 vertical levels (Neale et al., 2013).The Community Land Model (Lawrence et al., 2011) operates on the same horizontal grid as CAM4.The sea ice model is the Community Ice Code, version 4 (Hunke and Lipscomb, 2008), and the dynamic land ice component is inactive.The ocean physical model is identical to the ocean component of the Community Climate System Model version 4 (CCSM4) (Danabasoglu et al., 2012), except that shortwave absorption in the ocean is computed using prognostic chlorophyll fields, rather than a fixed satellitederived monthly climatology as in CCSM4.The ocean model has nominal 1 • horizontal resolution and 60 vertical levels.Mesoscale eddy transport is parameterized with an updated version of (Gent and McWilliams, 1990), where the eddyinduced advection coefficient, κ, is diagnosed as a function of space and time.Diapycnal mixing is represented using the K-Profile Parameterization of Large et al. (1994), and mixedlayer restratification by submesoscale eddies is parameterized using the method of Fox-Kemper et al. (2011).The biogeochemical-ecosystem ocean model incorporates multinutrient co-limitation on phytoplankton growth and specific phytoplankton functional groups (Moore et al., 2004(Moore et al., , 2013)), full carbonate system thermodynamics, sea-air CO 2 fluxes, and a dynamic iron cycle (Doney et al., 2006;Moore and Braucher, 2008).Phytoplankton calcification in the model is unaffected by variations in the saturation state of calcite or aragonite.Previous studies conducted with hindcast simulations of this model configuration reveal that the ocean physical state and air-sea CO 2 fluxes compare favorably with observations (Danabasoglu et al., 2012;Long et al., 2013).
Biogeochemical fields were initialized using data-based climatologies; for instance, DIC was from the Global Ocean Data Analysis Project (GLODAP; Key et al., 2004) and nutrients were from the World Ocean Atlas (Garcia et al., 2010).Subsequently, the fully coupled model was integrated for a period of 1000 years to allow the deep ocean to approach equilibrium; the tracer fields resulting from this spin-up procedure were used to initialize a 1000-year control simulation (Lindsay et al., 2014), in which atmospheric CO 2 was held constant at preindustrial levels (pCO atm 2 = 284.7 ppm).
By prescribing atmospheric CO 2 , this control simulation generates sea-air CO 2 flux variance that differs slightly from a control simulation using prognostic atmospheric CO 2 (Lindsay et al., 2014), owing to a lack of communication between land and ocean carbon reservoirs.Following Doney et al. (2006), variability in the control simulation is generated entirely from internal processes.During the first 100 years of the simulation, ocean [CO 2− 3 ] was not saved to disk, so our analysis is limited to the final 900 years of the simulation.Over this time, the global ocean drift in surface [CO 2− 3 ] is small (0.0029 mmol m −3 yr −1 ).

Model evaluation
Confidence in our interpretation of model output relies on the ability of the model to reproduce realistic estimates of mean surface [CO 2− 3 ] and its variability.We compare the annualmean surface [CO 2− 3 ] in the pre-industrial control simulation with reconstructed, pre-industrial surface [CO 2− 3 ] in Fig. 1.Pre-industrial surface [CO 2− 3 ] is reconstructed from observations of pre-industrial surface DIC (present-day DIC minus anthropogenic DIC) and present-day Alk as estimated by GLODAP (Key et al., 2004), combined with presentday estimates of temperature, salinity, silicate, and phosphate from the World Ocean Atlas (Locarnini et al., 2010;Antonov et al., 2010;Garcia et al., 2010) using a program developed for CO 2 system calculations (CO2SYS) with the preferred dissociation constants.Computing [CO 2− 3 ] this way is not strictly correct, owing to non-linear relationships between [CO 2− 3 ] and state variables, but no synthesized climatology of [CO 2− 3 ] from bottle sites exists.The model captures the large-scale distribution of surface [CO 2− 3 ] as estimated by the pre-industrial reconstruction (Fig. 1).Pre-industrial surface [CO 2− 3 ] is elevated in the tropical oceans, with the exception of the tropical Pacific cold tongue region, where persistent upwelling of DIC maintains low surface [CO 2− 3 ].High latitude, pre-industrial surface [CO 2− 3 ] is substantially lower than that in the tropics, due to upwelling of waters enriched in DIC and enhanced CO 2 solubility in these cold regions (Feely et al., 2004;Fabry et al., 2009).The model-estimated, pre-industrial surface [CO 2− 3 ] is noticeably lower than observed.The globally averaged surface [CO 2− 3 ] bias in the model is −16 mmol m −3 , excluding the Arctic Ocean and marginal seas.This model bias is caused by low biases in both surface Alk and DIC.The former bias is larger in magnitude and results from a combination of the prescribed carbonate dissolution profile, the representation of calcification, and a lack of riverine inputs in the model (Long et al., 2013).
We investigate the vertical distribution of salinitynormalized DIC (sDIC) and salinity-normalized Alk (sAlk) in the eastern Equatorial Pacific and compare with the vertical distribution of reconstructed pre-industrial sDIC and sAlk from observations in Fig. 2. As we will see later, the response of these tracers to ENSO variability depends on the their vertical gradients in this region.We define the eastern equatorial Pacific region according to the corresponding location of the biogeographical biome presented in Fay and McKinley (2014).In our model, the eastern equatorial Pacific is characterized by a large vertical gradient in sDIC, and a relatively small vertical gradient in sAlk, in general agreement with the observations in this region.Notably, the vertical gradient of modeled sDIC in the upper 1000 m is ∼ 25 % larger than that estimated from observations.Thus, the influence of ENSO on surface [CO 2− 3 ] is likely to be slightly exaggerated in the model.
We use the approach outlined in Friedrich et al. ( 2012) to evaluate simulated interannual variability and to compare it with observations collected at two open-ocean time series stations (BATS and HOT).At the model grid cells corresponding to the observational records, we generate probabil-ity density distributions of the standard deviation of annualmean surface [CO 2− 3 ] over a period of record sampled from the model that is the same length as the observational time series.At BATS, we calculate the standard deviation for 30 independent 30-year periods (Fig. 3a), and at HOT, we calculate the standard deviation for 36 independent 25-year periods (Fig. 3b).These are compared to the interannual standard deviation in de-trended, annual-mean surface [CO 2− 3 ] from the observational record, which is 30-years long at BATS and 25 years long at HOT (Bates et al., 2014).Figure 3 illustrates that the model underestimates interannual variability as compared to the observational data.This underestimate may be due to the mismatch in the spatiotemporal scales of model and observations; we would expect higher variance in the point-source observational record than in the smoothed model.As a result of the model's low variance bias, we expect our estimate of the detection timescale for the anthropogenic trend in [CO 2− 3 ] to be biased low, discussed next.

Surface carbonate variability and trend detection
We find the highest simulated interannual variability in preindustrial surface [CO 2− 3 ] in the eastern equatorial Pacific and at the boundaries between the subtropical and subpolar gyres (Fig. 4a).In the eastern equatorial Pacific, our model simulates a large vertical gradient in [CO 2− 3 ] (≈ Alk − DIC; Fig. 2), so variability is likely driven by changes in vertical motion here.The boundaries between the gyres exhibit the largest horizontal gradients in surface [CO 2− 3 ] (Fig. 1b), so changes in the shape of these gyres could have a large impact on local variability here (Friedrich et al., 2012).In the North Atlantic and North Pacific, the inter-gyre regions exhibit high variance on decadal (Fig. 4b) and longer timescales, as well, suggesting an influence from low-frequency climate variability.Finally, we note that simulated surface [CO 2− 3 ] exhibits low interannual variability in the Southern Ocean, consistent with previous modeling studies (Orr et al., 2005;Friedrich et al., 2012;Conrad and Lovenduski, 2015).How does the magnitude of the variability in surface [CO 2− 3 ] affect the detection of anthropogenic trends in aragonite ?We use monthly model output to investigate the length of the time series needed to detect an anthropogenic trend in aragonite with 90 % confidence, using the method of Weatherhead et al. (1998), At each location, we calculate the standard deviation (σ N ) and autocorrelation (φ) in the de-seasonalized, monthly anomalies of surface [CO 2− 3 ] and solve for the detection time (n * , Fig. 5a).This statistical technique has been applied successfully to ocean biogeochemical data in several previous studies (Henson et al., 2010;Beaulieu et al., 2013;Majkut et al., 2014;Lovenduski et al., 2015); it provides a way to quantify the importance of the variance and autocorrelation on the detection of the trend.It differs in practice from other statistical methods that estimate detectability or time of emergence of the trend (Ilyina et al., 2009;Friedrich et al., 2012;Hauri et al., 2013;Keller et al., 2014), in that it also includes the influence of temporal autocorrelation which affects sample size (Bretherton et al., 1999) and therefore trend detectability (Beaulieu et al., 2013).We note that the Weatherhead et al. ( 1998) method assumes a first-order auto-regressive process (AR(1)), and previous literature suggests that ocean carbon fluxes may be better represented by a higher-order auto-regressive process (e.g., AR(3); Séférian et al., 2013).As such, our detection times represent minimum values.
The estimate of the detection time is strongly influenced by the size of the anthropogenic trend to be detected (ω 0 ).Since n * is proportional to (1/ω 0 ) 2/3 , it takes longer to detect smaller trends and vice versa.As the simulation we analyze here is a pre-industrial control, carbonate chemistry is not influenced by anthropogenic factors, and we must turn elsewhere for an estimate of ω 0 .Bates et al. ( 2014 where [CO 2− 3 ] sat, aragonite is the carbonate ion concentration in equilibrium with mineral aragonite, which is primarily a function of pressure.Since variability in surface [CO 2− 3 ] sat, aragonite is an order of magnitude smaller than variability in surface [CO 2− 3 ] (not shown), we approximate the local, anthropogenic trend in surface [CO 2− 3 ] (ω 0 ) as the product of the global-mean, anthropogenic trend in aragonite (−0.0078 yr −1 ) and the local, model-estimated value of [CO 2− 3 ] sat, aragonite .The length of the time series needed to detect an anthropogenic aragonite trend of −0.0078 yr −1 (n * ) is shown in Fig. 5a.This detection time is spatially heterogeneous and ranges from 5 years at some locations to > 50 years at others.Equation (1) reveals that detection time is influenced by the ratio of the standard deviation to the trend (i.e., the noise-to-signal ratio, σ N /ω 0 , Fig. 6a) and the autocorrelation (φ, Fig. 6b).Further, the noise-to-signal ratio spatial pattern (Fig. 6a) is dominated by the noise (i.e., the variance, Table 1.The standard deviation (σ N ) and autocorrelation (φ) in the de-seasonalized monthly anomalies of surface [CO 2− 3 ] (mmol m −3 ) and pH, and the length of the time series in years (n * ) needed to detect an aragonite trend of −0.0078 yr −1 , and a pH trend of −0.0018 yr −1 with 90 % confidence at the stations discussed in Bates et al. (2014).Boldfaced values of n * indicate that the detection timescale is shorter than the length of the observational record at that location.Fig. 4a).Generally, we find long detection times in regions with high variance and high autocorrelation, such as at the boundaries between the subtropical and subpolar gyres in the Northern Hemisphere and in the Pacific cold tongue region, and short detection times in regions with low variance and low autocorrelation, such as in the Southern Ocean (cf. Figs. 6a and b, 5a).We note an additional area with high detection time in the eastern North Pacific Subtropical Gyre, adjacent to the California Current System, where variance is moderate, but autocorrelation is at a maximum.
[CO 2− 3 ] and aragonite are not directly measured in seawater, but rather derived from the measurement of other carbonate system variables, such as f CO 2 , pH, DIC, and Alk.It is thus of interest to know how the detection timescale for trends in [CO 2− 3 ] or aragonite compares to the detection timescale for trends in the measured parameters.Figure 5b shows the model-estimated detection time for a spatially uniform pH trend of −0.0018 yr −1 (the average pH trend of the seven time series analyzed in Bates et al. (2014), and the global-mean pH trend reported in Lauvset et al., 2015).As for aragonite , the pH trend detection time is elevated in the equatorial Pacific and inter-gyre regions, but we find much shorter detection times overall for the pH trend (cf.Fig. 5a and b). Figure 6 reveals that the shorter detection time for pH results from lower noise-to-signal ratios and lower autocorrelation than for [CO 2− 3 ].Thus, results from our model suggest that the anthropogenic trend in pH is detectable sooner than the anthropogenic trend in aragonite .These results imply that significant trends in pH will emerge before significant trends in aragonite at nearly every location in the ocean (with perhaps the exception of the Arctic Ocean).
We detail the model-estimated detection times for trends in aragonite and pH at the seven open ocean time series sites (Bates et al., 2014) in Table 1.For reference, the locations of the time series sites are shown as white circled xs in Fig. 5.The time series sites are ideally located in places with low variability in [CO 2− 3 ] and pH, so significant trends emerge more quickly there relative to other locations with high variance in these properties, such as the equatorial Pacific and inter-gyre regions (Fig. 4).According to our model calculations, the observational record is of sufficient length to detect the anthropogenic trends in aragonite and pH at all stations except Munida in the Southern Hemisphere, where the aragonite trend is not yet detectable (Table 1).Bates et al. (2014) report significant (at the 99 % level) decreases in aragonite and pH at all stations, with the exception of the Iceland Sea, where the trend in aragonite is not significant.That the model and observational records both suggest detectable anthropogenic trends at the time series sites is encouraging, but the agreement must be interpreted with some caution.We previously demonstrated that the model underestimates interannual variability at two of these sites, and our calculation of detection time from model output (see Eq. 1) assumes perfect temporal coverage at the sites.Taken together, the model detection time at the time series sites is likely to be an underestimate.

Modes and drivers of variability in surface carbonate
Given the important role of variance in trend detection, and the high magnitude of surface [CO 2− 3 ] variance in the Tropical Pacific, North Pacific, and North Atlantic (Fig. 4), it behooves us to characterize and understand the large-scale drivers of variability in these regions.Figure 7a shows the leading empirical orthogonal function (EOF) in tropical Pacific (18 • S-18 • N) surface [CO 2−  3 ].This EOF captures 55 % of the interannual variance in tropical Pacific surface [CO 2− 3 ] and its spatial pattern resembles that of ENSO SST variability.The wavelet power spectrum of the leading principal component (PC) in this region shows statistically significant spectral power on timescales associated with ENSO variability (3-7 years; Fig. 8a).We construct an ENSO index, defined as the annual-mean SST in a box bounded by 5 • S to 5 • N and 120-170 • W (Nino 3.4 region), from preindustrial control simulation model output.We note that, while the main characteristics of ENSO are well-captured by the model, the overall magnitude of ENSO in CCSM4 is overestimated (Deser et al., 2012b).Nevertheless, the modeled ENSO index is highly correlated with the leading PC of surface [CO 2− 3 ] (r = −0.94).This suggests that ENSO exerts a strong control on surface [CO 2− 3 ] variability in the tropical Pacific.
We investigate the drivers of variability in tropical Pacific surface [CO (S) using a linear Taylor expansion, where the terms are the perturbations in state variables associated with ENSO, which we estimate using the regression coefficients of the local value with the leading PC of surface [CO 2− 3 ] in the tropical Pacific region, and the partial derivatives are estimated from CO2SYS using finite difference approximation.We separate the contribution from freshwater (fw) fluxes on DIC and Alk by expanding the derivatives to include sDIC and sAlk (see Appendix for more details).Sub- whose individual terms are shown in Fig. 7. Crosscorrelations among the variables in the Taylor sum may lead to imprecise results in some locations.
Our analysis demonstrates that variations in sDIC are the primary driver of [CO 2− 3 ] variability in the tropical Pacific, with smaller, opposing contributions from sAlk and freshwater (Fig. 7).Contributions from temperature and salinity play comparatively smaller roles.As the east/central equatorial Pacific is associated with a large vertical gradient in sDIC (Fig. 2), variations in vertical exchange here likely dominate surface sDIC variability.El Niño events are associated with relaxed trade winds and less upwelling of DIC-rich water in the east/central equatorial Pacific, raising the surface [CO 2− 3 ]; the opposite is true during La Niña events, when we would expect anomalously low surface [CO 2− 3 ] in these regions.In our simulation, the ENSO index has a strong negative (positive) correlation with surface sDIC ([CO 2− 3 ]) in the east/central equatorial Pacific (not shown).
The leading EOF of surface [CO 2− 3 ] variability in the North Pacific (20-70 • N) is characterized by a broad, horseshoe-shaped area of [CO 2− 3 ] off the coast of North America that is out of phase with [CO 2− 3 ] in the western subtropical gyre (Fig. 9a).This EOF explains 35 % of the variance in surface [CO 2− 3 ] in the North Pacific basin and its spa-tial pattern resembles the PDO.The wavelet spectrum for the associated first PC (Fig. 8b) reveals high spectral power at 20-30 year periods, similar to the PDO, but also detects significant spectral power at high frequencies (3-7 years), and low frequencies (> 60 years).We define a model PDO as the leading EOF of sea surface temperature anomalies in the North Pacific (20-70 • N); it has a generally realistic spatial pattern and amplitude (Deser et al., 2012b).The model PDO is highly correlated with the leading PC of surface [CO 2− 3 ] in the North Pacific (r = −0.94),suggesting that the PDO has a large influence on surface [CO 2− 3 ] in this region.We decompose the leading EOF of surface [CO 2− 3 ] in the North Pacific into contributions from sDIC, sAlk, freshwater, temperature, and salinity using Eq. ( 4).Our analysis reveals that sDIC contributes to most of the variability in surface [CO 2−  3 ] in this region, with smaller, opposing contributions from sAlk and freshwater, and near-zero contributions from temperature and salinity (Fig. 9).Positive/warm phases of the PDO are associated with southerly winds off the coast of California, suppressing the upwelling of carbon-rich water and elevating surface [CO 2− 3 ], while the western subtropical gyre experiences increased upwelling of DIC and lower surface [CO 2− 3 ].The opposite is true during negative/cold phases of the PDO, when we find anomalously high (low) surface [CO 2− 3 ] in the eastern (western) North Pacific.We find a significant positive correlation between the PDO and surface [CO 2− 3 ] off the coast of California and a significant negative correlation between the PDO and surface [CO 2− 3 ] in the western North Pacific subtropical gyre (not shown).
In the North Atlantic (40-80 • N), the leading EOF of surface [CO 2− 3 ] explains 18 % of the variance and is characterized by a tripole pattern, with nodes centered in the subpolar gyre, the inter-gyre region, and the subtropical gyre (Fig. 10).The wavelet time-frequency diagram of the lead- ing PC (Fig. 8c) looks like a red noise process (Torrence and Compo, 1998) with significant spectral power at 20-40 year and longer periods.A Taylor series decomposition of the variability in this region using Eq. ( 4) reveals that, similar to the other two regions, sDIC variability is the dominant driver of [CO 2− 3 ] variations here, with smaller, opposing contributions from sAlk, freshwater, and salinity in the inter-gyre region (Fig. 10) where variance in [CO 2− 3 ] is at a maximum (Fig. 4).The leading PC of [CO 2− 3 ] in this region is weakly correlated with the NAO (r = −0.16)and the Atlantic Multidecadal Oscillation (AMO, r = −0.33).Positive phases of the NAO and AMO are associated with anomalously high [CO 2− 3 ] in the subpolar and subtropical gyres and anomalously low [CO 2− 3 ] in the inter-gyre region, though the weak correlations preclude further investigation of the mechanistic link between these modes of climate variability and surface [CO 2− 3 ] in this region.

Conclusions
We analyze output from a 1000-year, pre-industrial control simulation of an Earth System Model in order to quantify, characterize, and understand the drivers of variability in the surface ocean carbonate ion concentration.We find the highest variability in the tropical Pacific, where the model simulates large vertical gradients in [CO 2− 3 ], and at the boundaries between the subtropical and subpolar gyres in the Northern Hemisphere.High variance coupled with high autocorrelation results in long detection times (> 50 years) for anthropogenic trends in aragonite at these locations, whereas we find short detection times (∼ 15 years) at the open ocean time series sites and in the Southern Ocean, where variance is at a minimum.Our results suggest that the detection timescale for anthropogenic trends in pH is shorter than that for aragonite , owing to smaller noise-to-signal ratios and lower autocorrelation in pH.We further characterize the surface [CO 2− 3 ] variance in the Tropical Pacific, North Pacific, and North Atlantic using EOF and wavelet analysis.The leading mode of variability in the tropical Pacific exhibits statistically significant spectral power at 3-7 year timescales and is highly correlated with ENSO.In the North Pacific, the leading EOF of [CO 2− 3 ] has the characteristic horseshoe pattern of the PDO, with peak spectral power at 20-30 year timescales.In the North Atlantic, the wavelet decomposition of the leading mode of variability appears similar to a red noise process, with weak correlations with the NAO and AMO.In all locations, surface [CO 2− 3 ] variability is driven by variability in sDIC, with smaller opposing contributions from sAlk and freshwater dilution.Temperature and salinity variability contributes very little to [CO 2− 3 ] variance.In all regions, climate variability imposes changes in ocean circulation that likely mediate the vertical and lateral advection of DIC into the surface ocean.
While the spatial pattern of surface [CO 2− 3 ] variability in our model is similar to that reported in other studies (e.g., Friedrich et al., 2012), a recent study of natural carbon uptake variability from centuries-long simulations of six ESMs suggests that we should expect the magnitude of this variability to differ from model to model (Resplandy et al., 2015).Thus, the detection times for trends in [CO 2− 3 ] and pH, which themselves are a function of the variability, are also likely to differ from model to model.Future investigations characterizing the uncertainty in internal variability from ESM ensembles are therefore needed.
Our results provide meaningful perspective on the trends in aragonite and pH reported in the literature.Our statistical analysis of model output suggests that, due to low variance in [CO 2− 3 ] and pH, significant anthropogenic trends in aragonite and pH are already or nearly detectable at the open-ocean time series sites, whereas high variability and autocorrelation in [CO 2− 3 ] may obscure the detection of anthropogenic trends in aragonite in the tropical Pacific, North Pacific, and North Atlantic.Our characterization of the spatial pattern and frequency of natural [CO 2− 3 ] variability in these high variance regions suggests that it is largely driven by large-scale modes of internal climate variability, such as ENSO, PDO, NAO, and AMO.One should consider the phasing of these modes of climate variability, therefore, when interpreting trends calculated from carbonate chemistry data collected in these regions.
Ongoing efforts aimed at optimizing the ocean acidification observational network (e.g., the Global Ocean Acidification Observing Network, GOA-ON) stand to benefit from the results presented here.Our analysis suggests that frequent, sustained observations are required to capture the natural variability that is key to detecting significant anthropogenic trends in surface aragonite and pH.In regions where variability is muted (e.g., the subtropical gyres), detectable trends can emerge from sparse or short records, whereas, in regions with heightened variability (e.g., the equatorial Pacific), the observational record must be denser and longer in order to detect the same trend.Further, our finding that the detection timescale for pH is shorter than that for aragonite nearly everywhere highlights the need for continued underway and float-based pH measurements.
et al.: Surface ocean carbonate variability tions.In the equatorial Pacific,

Figure 3 .
Figure 3. Blue bars: probability density of interannual standard deviation of simulated [CO 2− 3 ] at (a) BATS, based on 30 independent 30-year periods and (b) HOT, based on 36 independent 25-year periods.Black line: interannual standard deviation of surface [CO 2− 3 ] at (a) BATS and (b) HOT, based on detrended annual-mean observations.

Figure 5 .
Figure 5. Model-estimated length of time series in years needed to detect (a) an aragonite trend of −0.0078 yr −1 , and (b) a pH trend of −0.0018 yr −1 with at least 90 % confidence.White circled xs indicate the locations of the time series stations discussed in Bates et al. (2014).

Figure 9 .
Figure 9.As in Fig. 7, but for the North Pacific.

Figure 10 .
Figure 10.As in Fig. 7, but for the North Atlantic.