Interactive comment on “ Climate impacts on multidecadal p CO 2 variability in the North Atlantic : 1948 – 2009

This manuscript attempts to establish a statistical and conceptual link between AMO, SST, ocean pCO2 and air-sea CO2 fluxes. The primary effects come from the SST and the vertical transport of DIC as explained by the previous publications from the same group. What is significant in this new paper would be further linking this to the AMO on multi-decadal timescales. The manuscript is written in a reasonably coherent and clear way. While I generally support the qualitative conclusions, I have several suggestions as listed below.


Introduction
To date, the ocean has removed approximately one-third of all anthropogenic carbon emitted to the atmosphere and has, thus, substantially damped climate warming (Khatiwala et al., 2009;Sabine et al., 2004).As carbon dioxide emissions continue to increase due to fossil fuel emissions and cement production, there is significant interest in better understanding the ocean carbon cycle.Due to the limited instrumental record and sparse data, multidecadal variability of the ocean carbon sink remains poorly constrained.The North Atlantic, in particular, is a region of highly concentrated carbon uptake (Takahashi et al., 2009) and of significant carbon cycle variability related to variations in the climate, with multiple studies finding an association with the North Atlantic Oscillation (NAO; Fay and McKinley, 2013;Schuster et al., 2013;Terry, 2012;Levine et al., 2011;McKinley et al., 2011;Löptien and Eden, 2010;Ullman et al., 2009;Thomas et al., 2008).However, data are sparse, processes are complex, and the timescales for studies have differed, and this has complicated a clear elucidation of the mechanisms of North Atlantic carbon cycle variations.Schuster et al. (2009) analyzed in situ partial pressure of CO 2 (pCO 2 ) measurements and suggested a substantial decline in North Atlantic carbon uptake from the mid-1990s to the mid-2000s.Le Quéré et al. (2010) also interpreted observations and models to conclude that there had been a decline in the North Atlantic sink from 1981 to 2007 due to changing wind patterns and increasing sea surface temperature (SST).Metzl et al. (2010) focused on subpolar surface ocean carbon cycle changes between 1993 and 2008, and also concluded that there had been a reduction in carbon uptake.In situ pCO 2 measurements have also been synthesized to illustrate the strong sensitivities of such changes to the locations and time frame chosen for the analyses (Fay and McKinley 2013;McKinley et al., 2011).The substantial spatial heterogeneity and temporal variability in the North Atlantic complicates efforts to use sparse observations to quantify carbon uptake.Thus, the magnitude and mechanisms affecting North Atlantic carbon cycle variability remain loosely constrained.The present study takes advantage of the full spatial and temporal coverage of a regional numerical model to gain new in-sights into the mechanisms of variability of North Atlantic pCO 2 .
As shown by Ullman et al. (2009) in a 15-year simulation (1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006), internal variability in the North Atlantic is partially obscured by the large, quasilinear trend of CO 2 flux into the ocean that is driven by increasing CO 2 emissions.To examine the carbon sink variability that is partially masked by this large carbon influx, we use a hindcast model from 1948 to 2009 forced with the preindustrial atmospheric CO 2 concentration and realistic climate.As described below, we find that the basin-average SST is associated with the leading mode of surface ocean pCO 2 variability.This SST signal, in turn, includes an upward trend due to greenhouse gas emissions and a signal of internal variability characterized by the Atlantic Multidecadal Oscillation (AMO; Kerr, 2000).

Physical-biogeochemical-ecosystem model
The MIT Ocean General Circulation Model (Marshall et al., 1997a, b) has been regionally configured for the North Atlantic between 20 • S and 81.5 • N (Bennington et al., 2009;Ullman et al., 2009).The model has a horizontal resolution of 0.5 • latitude and 0.5 • longitude and 23 vertical levels beginning with a resolution of 10 m thickness at the surface and increasing to 500 m thickness at depths greater than 2200 m.The Gent-McWilliams (Gent and McWilliams, 1990) eddy parameterization and the K-profile parameterization (KPP) boundary layer mixing scheme (Large et al., 1994) were employed to model sub-grid-scale processes.Daily fields from NCEP/NCAR Reanalysis I force the model from 1948 to 2009 (Kalnay et al., 1996).SST and SSS (sea surface salinity) are relaxed to monthly historical SST (Had1SSTv1.0,Rayner et al., 2003) and climatological SSS (Antonov et al., 2006) observations, on the timescale of 2 and 4 weeks, respectively.Glacier melt and/or river discharge are not included in the model forcing; instead the SSS relaxation approximates these impacts.Freshwater (evaporation−precipitation) forcing and SSS relaxation impact both salinity and tracer concentrations.In lieu of an active sea ice simulation, observed fractional ice from NCEP Reanalysis 1 is applied with interpolation to daily resolution.
For open boundary conditions, a sponge layer exists at 20 • S, and over the first 5 • of latitude to the north there is restoration to climatological temperature, salinity, dissolved inorganic carbon (DIC), and phosphate fields.For temperature and salinity, there is also a sponge layer at Gibraltar.More discussion of the sponge layer can be found in Ullman et al. (2009).
Carbon (inorganic and dissolved, and particulate organic), alkalinity (ALK), phosphorus, silica, and iron cycling are explicitly included in the biogeochemical model.Carbonate chemistry is modeled as in Follows et al. (2006).The objective of this simulation is to identify climate impacts on the natural carbon cycle without the complication of the large CO 2 flux into the ocean that is observed.Thus, atmospheric pCO 2 is fixed at a constant, preindustrial level of 278 ppmv.
The physical model was spun up for 100 years.The observed meteorological forcing for 1948-1978 was cycled three times (90 years), and then 1948 was repeated 10 years prior to beginning the simulation analyzed here.Following the physical spin-up, the biogeochemical model was initialized using preindustrial estimates for DIC and ALK climatology from the Global Data Analysis Project (GLODAP) database (Key et al., 2004).The biogeochemical model was then spun up for an additional 100 years, long enough to eliminate drift in the biogeochemical parameters.The percent change over the last 5 years of spin-up in the basinaveraged surface DIC field is 0.00046 % per year.For comparison, the percent change in DIC from a high AMO (1955) to low AMO (1975) is 0.012 % per year, two orders of magnitude greater than drift at the end of the spin-up.This indicates that a 100-year biogeochemical spin-up is sufficient to eliminate model drift that would impact our upper-ocean analysis.Following spin-up, the model was then run with NCEP/NCAR daily forcing fields for 1948-2009.
Model physics across the North Atlantic as well as pCO 2 , DIC, and ALK at the Bermuda Atlantic Time Series (Bates, 2007) and in the subpolar North Atlantic have been compared to results from a previous simulation using the same model forced with observed atmospheric pCO 2 for 1992-2006 (Ullman et al., 2009).Comparison of this simulation to estimates of the preindustrial vertical profile of DIC in the subpolar gyre indicates good performance by the model (Supplement Fig. S1 and text).Mikaloff-Fletcher et al. (2007) estimated the preindustrial, or "natural", air-to-sea CO 2 flux in the North Atlantic with an ocean inversion that incorporated climatological circulations estimated from 10 ocean circulation models.For the North Atlantic from 0 to 75 • N, they find an uptake of 0.27 ± 0.07 PgC yr −1 .The mean natural CO 2 flux averaged over the same spatial domain in our simulation is consistent, 0.23 PgC yr −1 .In total, our comparison to available data indicates that the model is capable of robustly simulating the carbon biogeochemistry of the North Atlantic and its response to climate variability.

Post-processing
CO 2 flux into the ocean is proportional to the partial pressure difference between the atmosphere and ocean surface: In this analysis, we can directly relate higher pCO ocn 2 to a reduction in CO 2 flux, since atmospheric pCO 2 is fixed.pCO 2 variability sets the sign and magnitude of flux changes on both seasonal and interan- nual timescales (Takahashi et al., 2009;Watson et al., 2009, Le Quéré et al., 2010).pCO 2 is decomposed into contributions from temperature and chemical effects using model output and the full carbonate equations (Follows et al., 2006).As in Ullman et al. (2009), pCO 2 -SST is found by allowing only SST to vary in the full carbonate equations for pCO 2 ; i.e., all other variables (DIC, ALK, SSS, phosphate, silica) are held constant at their long-term mean values.pCO 2 -chem is found by holding SST constant and allowing the rest of the input variables to vary; for pCO 2 -DIC, only DIC varies.
Model diagnostics for DIC are the monthly mean tendency terms (in mmol m −3 yr −1 ) due to individual modeled processes and are calculated at each time step during the model simulation (Ullman et al., 2009).Monthly mean diagnostics for the surface layer DIC change due to horizontal and vertical advection and diffusion, net biological processes (pri-mary production and respiration), freshwater input/removal, and air-sea CO 2 flux are used.
The AMO index for the model is calculated using modeled SST and observed global Had1SSTv1.0(Rayner et al., 2003) using the approach of Wang and Dong (2010).This approach regresses the area-weighted global mean Had1SST time series onto area-weighted basin-wide mean North Atlantic SST time series (NASST).This regressed index is subtracted from the total NASST to define the AMO.The combined SST signal is, thus, decomposed into contributions from globally increasing SST (SST trend) and the internal variability of the AMO (Fig. 1d).In order to focus on the decadal timescale variability, all time series are smoothed with a standardized 121-month box smoother.

Multidecadal variability
To determine the leading mode of variability in surface ocean pCO 2 , principle component analysis is employed.The first empirical orthogonal function (EOF1) patterns and smoothed principle components (PCs) for monthly, 13month smoothed total pCO 2 , and the SST contribution to pCO 2 (pCO 2 -SST) are shown in Figure 1a-c.To determine the change in pCO 2 anomalies described by EOF1 at a specific point in time, the value of the PC1 at that time can be multiplied by the EOF1 pattern.The percent of variance in the total field explained by the EOF1 pattern is 18 and 38 % for pCO 2 and pCO 2 -SST, respectively.In both cases, the EOF1 patterns are statistically distinct from their EOF2 patterns, which are discussed in Sect. 4. This EOF analysis unveils the basin-scale coherent variability.There is remaining variability in coherent secondary large-scale modes (e.g., EOF2) or at scales smaller than the whole basin.That large-scale modes of climatic variability tend to capture 10-40 % of variance has been documented across many climate variables, including global SST and tropospheric winds (von Storch and Zwiers, 1999), Southern Ocean geopotential heights (Thompson and Wallace, 2000), and pCO 2 throughout the Pacific (McKinley et al., 2004(McKinley et al., , 2006)).That EOF1 of pCO 2 captures the patterns of multidecadal large-scale change is further evidenced by plots of 20-year anomalies of pCO 2 (Fig. S2).
The correlation between PC1-pCO 2 and the areaweighted, basin-averaged SST is 0.88 (Fig. 1c, Table S1 in the Supplement).An increase in temperature increases pCO 2 by reducing solubility, which is illustrated by the pCO 2 -SST EOF1 pattern.PC1-pCO 2 and PC1-pCO 2 -SST are highly correlated (Fig. 1c, r = 0.91) but have distinct EOF1 patterns, particularly in the subpolar gyre (Fig. 1a, b).This is consistent with the pCO 2 in the subpolar gyre also being significantly impacted by changes in DIC supply which, in turn, are associated with the AMO.EOF1 for pCO 2 -chem  and pCO 2 -DIC explain 32 and 25 % of the variance, respectively (Fig. 2a, b), and these PC1s are highly correlated with the AMO: r = 0.99 and 0.96, respectively (Fig. 2d, Table S1).
Alkalinity can also affect pCO 2 -chem since increased alkalinity reduces pCO 2 .PC1 for EOF1 of pCO 2 -ALK (Fig. 2c) does not correlate highly to PC1s of total pCO 2 or pCO 2 -chem (r = −0.25;r = 0.44, respectively), or to the AMO (see Table S1).Though alkalinity does contribute to the spatial pattern shown in the EOF1 of pCO 2 -chem, the temporal evolution of this pattern differs substantially and is not strongly connected to the AMO or to EOF1 of pCO 2 .Therefore, we focus on the more direct relationship between pCO 2 -DIC and pCO 2 -chem for the rest of the paper, and reserve the alkalinity relationships for future in-depth analysis.
The AMO, an index of internal North Atlantic SST variability, declines (cools) until 1975 and rises thereafter (Fig. 1d).Taking the last half of the time series as an example, increasingly positive AMO corresponds to a decrease in pCO 2 -chem, with the strongest declines in the subpolar gyre and driven by reduced pCO 2 -DIC (Fig. 2).This occurs in opposition to the direct effect on pCO 2 of warmer NASST (Fig. 1b, c), driven jointly by the increasingly positive AMO and the warming trend (Fig. 1d).SST and chemical terms vary inversely because higher SST (AMO+) enhances stratification, leading to a shoaling of mixed-layer depths (MLDs) over most of the gyre (discussed in Sect.4).This shoaling in turn limits the amount of deep, carbon-rich water that is mixed to the surface, reducing pCO 2 -DIC and pCO 2 -chem (Ullman et al., 2009).The correlation of PC1-pCO 2 -chem and PC1-pCO 2 -DIC with PC1-pCO 2 -SST is 0.90 and 0.91, respectively (Table S1).Mechanisms of AMO impacts on pCO 2 -chem in the subpolar gyre will be explored further below.

Regression analysis
Regression of the AMO, SST trend, and total SST (Fig. 1d) onto monthly pCO 2 , pCO 2 -SST, and pCO 2 -chem further illustrates that temperature and chemical responses tend to act in opposition to one another, damping total pCO 2 responses across the basin (Fig. 3).This analysis complements the above EOF analysis by allowing the use of the same index of temporal variability across all fields.Previous studies with observations and models have shown that pCO 2 -chem dominates the seasonality of pCO 2 in the subpolar gyre, via strong vertical supply of DIC in winter that drives up pCO 2 and biological DIC drawdown in summer that reduces pCO 2 .Temperature impacts oppose these seasonal oscillations but are of weaker amplitude (Körtzinger et al., 2008;Takahashi et al., 2002).Models have shown similar opposing influences with respect to interannual variability (Ullman et al., 2009;McKinley et al., 2004).These regressions illustrate that positive AMO leads to higher pCO 2 -SST throughout the basin (Fig. 3b).The response is strongest north of 35 • N with a clear maximum to the east of Newfoundland.Simultaneously, positive AMO is associated with a reduction in pCO 2 -chem (Fig. 3c).The pCO 2 -chem signal is also strongest to the north.The overall effect is a decrease in total pCO 2 north of 45 • N and a slight increase in pCO 2 in the eastern subtropical gyre (Fig. 3a).
When responding to the global SST trend, pCO 2 -SST more heavily controls the response of the total pCO 2 field (Fig. 3d, e).The pCO 2 -SST response is strongest along the Gulf Stream and east of Newfoundland, and it also increases somewhat off the coast of Europe and Africa.pCO 2 -chem exhibits some decline in the Gulf Stream region and has a small response elsewhere (Fig. 3f).
Regression with the total NASST time series (Fig. 1) illustrates the combined effects of the AMO and trend signals (Fig. 3g-i).A positive anomaly of NASST depresses total pCO 2 in the subpolar gyre, consistent with the AMO impact found above.Positive NASST also increases total pCO 2 off north Africa, consistent with the impact of the SST trend.pCO 2 -SST both increases off Africa and has a strong maximum in the Gulf Stream region east of Newfoundland with positive NASST anomalies.The pCO 2 -chem response is slightly weaker in the subpolar gyre than for the AMO alone (Fig. 3a, i).

DIC diagnostics
To further investigate the chemical term response to the AMO, model diagnostics for the DIC field are regressed upon the AMO index.Diagnostics are modeled rates of change in DIC due to one of five processes that have been saved at every model time step.Physical processes are separated into horizontal advection and diffusion (DIC-horz) and vertical advection and diffusion (DIC-vert).DIC-phys is the sum of vertical and horizontal transport, showing the net effect of physical transport on DIC (Fig. 4).The rate of DIC supply is also affected by biological processes involving DIC incorporation into organic matter and remineralization back to inorganic (DIC-bio), net precipitation/evaporation that dilutes or concentrates DIC (DIC-fresh), and the air-sea flux of CO 2 (DIC-flx) (Fig. 5).The focus on DIC is justified by the fact that pCO 2 -chem change has the same pattern and is highly correlated with pCO 2 -DIC change (Fig. 2).The focus on the AMO is justified by its strong imprint on pCO 2 through pCO 2 -chem (Figs. 2, 3).
For the long-term average, vertical advection and diffusion are positive along the Gulf Stream and in the subpolar gyre due to deep winter MLDs that mix up high-DIC water www.biogeosciences.net/13/3387/2016/Biogeosciences, 13, 3387-3396, 2016 from below (Fig. 4a).Horizontal DIC advection and mixing removes this vertically supplied DIC along the Gulf Stream and in the western subpolar gyre (Fig. 4b).While the vertical and horizontal components tend to have opposing influences, the net effect is a positive DIC supply to the subpolar gyre, as shown by mean DIC-phys (Fig. 4c).With positive AMO, vertical advective and diffusive fluxes of DIC decrease in the Irminger Sea and Iceland Basin, while they increase in the Labrador Sea and east of Newfoundland (Fig. 4d).These changes are consistent with AMO-related MLD changes outside of the Labrador Sea (Fig. 6) and changes in the basinscale barotropic stream function indicating a weakened subpolar gyre (Fig. 7).The effect of this is to shift the central DIC-vert maximum to the west.With positive AMO, horizontal advection and diffusion largely respond to changes in vertical advection and diffusion, with less horizontal divergence (a positive change) in regions where the vertical supply is reduced (Fig. 4e).The net effect shown by DIC-phys reveals an overall reduction in DIC supply (Fig. 4f), consistent with a weaker subpolar gyre circulation and shallower MLDs that reduce the vertical supply of DIC.Hakkinen and Rhines ( 2009) illustrate and increased penetration of subtropical waters into the subpolar region from the 1990s to the 2000s, consistent with a weaker subpolar gyre circulation.
The changes in MLD and stream function are also in agreement with results from Zhang (2008), who links the observed spin-down of the subpolar gyre in the 1990s to an enhanced Meridional Overturning Circulation (MOC), using a combination of satellite altimeter observations and results from a 1000-year coupled ocean-atmosphere model simulation.
Mean DIC impacts from physics, biological processes, freshwater, and air-sea flux are shown in Fig. 5a-e.The net impact of biology is to remove DIC from the surface of most of the region, with the most intense removal along the Gulf Stream (Fig. 5a).The smaller impact of evaporation and precipitation is to concentrate DIC in the subtropics and to dilute it in the subpolar gyre (Fig. 5b).The air-sea CO 2 flux term is also small, positive north of about 35 • N and negative to the south (Fig. 5c).AMO-related change in the biological removal of DIC indicates additional removal (negative anomaly) occurring in the same region where horizontal flux increases, consistent with biological stimulation through an increased supply of nutrients from the subtropical subsurface along the "nutrient stream" (Williams et al., 2006).There is reduced biological productivity, and thus a reduction of DIC loss (a positive DIC anomaly), in other parts of the basin, which is consistent with satellite observations from the late 1990s to the mid-2000s (Behrenfeld et al., 2006).Changes in surface ocean DIC content due to freshwater fluxes and air-sea CO 2 flux with the AMO are small.Across the basin, the net DIC change associated with AMO is negative, with the strongest negative changes occurring in the subpolar gyre (Figs.2b, 3c)

Discussion and conclusions
In this North Atlantic regional model forced with preindustrial pCO 2 and realistic climate from 1948 to 2009, SST is the dominant driver of pCO 2 variability, with both longterm anthropogenic warming and the AMO playing important roles.The AMO strongly influences chemical change, which in turn is mostly driven by DIC.DIC changes, in turn, are due primarily to changes in vertical and horizontal advecwww.biogeosciences.net/13/3387/2016/Biogeosciences, 13, 3387-3396, 2016 tion and mixing.Changing biology has the most important secondary effect and largely damps the anomalies caused by advection and mixing.Freshwater and CO 2 fluxes changes are slight.
Our findings linking the AMO to natural carbon cycle variability in the North Atlantic are consistent with the study of Séférian et al. (2013), who also found an AMO-like signal dominated North Atlantic pCO 2 variability in a 1000-year Earth system model simulation with constant pCO 2 .Other studies have focused on the relationship between the North Atlantic Oscillation and CO 2 flux using models and observations (Löptien and Eden, 2010;Ullman et al., 2009;Schuster et al., 2009;Thomas et al., 2008).Consistent with these previous studies, the NAO is the second mode of variability in this simulation (Figs.S4, S5), and the corresponding principle components are highly correlated with the NAO (Table S2).The shorter time frame for most previous studies explains, in part, the difference in attribution to the AMO as opposed to NAO.Our results are broadly consistent with previous studies in the finding that physical variability is the dominant driver of variability in the North Atlantic surface ocean carbon cycle.
The NAO and AMO may, in fact, be linked through the MOC, with a positive NAO enhancing MOC, which over time warms SSTs and leads to a positive AMO.Stronger overturning (enhanced MOC) may result in positive SST anomalies and thus a positive AMO phase over time, with some lag between the peak MOC and peak AMO response.Latif et al. (2004) used a model to show that oceanic heat transport related to the MOC leads the SST response (potentially the AMO) by about 10 years, while Delworth and Mann (2000) found a 10-year lag between subsurface temperature anomalies and MOC but less of a lag with SST.The precise mechanisms remain in debate due to different model findings and a lack of observational constraints (Delworth and Mann, 2000;Knight et al., 2005;Dima and Lohmann, 2007;Latif et al., 2006;Kavvada et al., 2013, and references therein).In our simulation, the NAO and MOC are significantly correlated (r = 0.57, Table S1) and there is also a high correlation (r = 0.86) between the NAO (Fig. S3) and the 15-year-lagged AMO.These correlations are broadly consistent with the above-postulated NAO-MOC-AMO relationship.On the other hand, Booth et al. (2012) suggest that the AMO may be driven, in fact, by atmospheric aerosol variability, so it is possible that there is no such AMO-MOC relationship at all.Future modeling and observations should further elucidate these connections, which reach beyond the scope of this study.
We find multidecadal variability in the natural carbon cycle of the surface North Atlantic to be dominated by an SST trend and multidecadal SST variations captured by the AMO index.Variability linked to the AMO influences both pCO 2 -SST and pCO 2 -chem.In the subpolar gyre, the positive SST influence on pCO 2 is overwhelmed by reduced supply of DIC to the surface ocean through mixing and advection, the net impact being reduced pCO 2 .The reduction in mixing is associated with shoaling of MLDs and a weaker subpolar gyre circulation, both associated with warmer SSTs (positive AMO).In the subtropics, the SST impact is stronger, and thus pCO 2 is increased under the influence of positive AMO and positive SST trend.
These findings are consistent with observed relationships between trends in surface ocean pCO 2 and trends in atmospheric pCO 2 since the 1980s (Fay and McKinley, 2013).In the North Atlantic subpolar gyre, trends in surface ocean pCO 2 lagged the trend in atmospheric pCO 2 from the early to mid 1990s to the late 2000s, which is consistent with the AMO and the SST trend reducing DIC supply to the subpolar gyre as found in this study.On smaller spatial scales and in shorter time frames, trends in ocean pCO 2 can differ (Fay and McKinley, 2013;Metzl, 2010;Watson et al., 2009, Schuster et al., 2009), which can be reasonably attributed to shorter-term and smaller spatial scale variability.We also find that warming has contributed to the observed pCO 2 increase from the 1980-90s through the 2000s throughout the basin.These model results allow a mechanistic attribution of these observed changes in North Atlantic pCO 2 to the com-bined effect of the AMO and a positive SST trend due to anthropogenic climate change.
The Supplement related to this article is available online at doi:10.5194/bg-13-3387-2016-supplement.
Figure 1.(a) EOF1 of total pCO 2 (µatm) and (b) EOF1 of pCO 2 -SST (µatm), explaining 18 and 38 % of total variance, respectively.(c) PC1-pCO 2 (orange); PC1-pCO 2 -SST (pink); and areaweighted, basin-averaged standardized North Atlantic SST time series (black).(d) Area-weighted, basin-averaged (0-70 • N, 98 • W-19.5 • E) North Atlantic SST from Had1SST (black); global areaweighted SST regressed onto North Atlantic SST (blue); and AMO index created by subtracting the global regression from the North Atlantic SST (red).All indices are standardized by 1σ .Time series smoothed with a 121-month box smoother.Two small coastal areas off Africa and South America were excluded in (a) and (b) due to the presence of localized, anomalously strong upwelling in the early 1960s that precluded elucidation of the large-scale pattern.