Dynamical and biogeochemical control on the decadal variability of ocean carbon fluxes

Introduction Conclusions References


Introduction
In recent years, several observation-based and model-based studies have suggested that ocean anthropogenic carbon uptake has slowed down due to the impact of anthropogenic forced climate change in various ocean regions (e.g.Corbière et al., 2007;Le Quéré et al., 2007).However, the fact that these changes in ocean carbon fluxes are caused by anthropogenic forced climate change within the last decades is still debated (e.g.Ullman et al., 2009;McKinley et al., 2011).
Among these regions potentially impacted by anthropogenic forced climate change, the North Atlantic has been extensively studied over the three last decades (Watson et al., 2011;Schuster et al., 2009;Brown et al., 2010;Metzl et al., 2010;Corbière et al., 2007;McKinley and Follows, 2004;McKinley et al., 2011).Debate prominently arises from disagreement between causes of trends in the North Atlantic carbon sink estimated from data or numerical model output.On the basis of in situ observations, several authors (i.e.Corbière et al., 2007;Watson et al., 2011;Brown et al., 2010;Metzl et al., 2010) have reported a more rapid growth in surface ocean pCO 2 than in atmospheric pCO 2 , which translates into a decline of the net ocean carbon uptake over the last decade.Possible mechanisms driving this decline are associated to climate-induced modifications of both oceanic and atmospheric dynamics, e.g., rising sea-surface temperature (Corbière et al., 2007), increasing ocean stratification (Schuster et al., 2009) or changes in horizontal oceanic currents owing to a shift in the North Atlantic Oscillation (Thomas et al., 2007;Schuster et al., 2009).However, model R. Séférian et al.: Decadal variability of ocean CO 2 fluxes reanalysis studies, based on GFDL, CCSM and MIT ocean general circulation reanalyses, have demonstrated that the weakening of the North Atlantic carbon sink over the last three decades (i.e., 1970s to 2000s) cannot be really attributed to anthropogenic climate change, but rather to the natural variability of the ocean carbon sink driven by the North Atlantic Oscillation (McKinley et al., 2011;Thomas et al., 2008;Löptien and Eden, 2010) as a part of the Northern Annular Mode (NAM).
Similarly to the North Atlantic, many aspects of the Southern Ocean dynamics and carbon cycle have exhibited trends over the last decades.Model simulations and inversions of atmospheric data have suggested that the southward shift of westerlies concomitantly with their strengthening had dampened the Southern Ocean carbon sink (e.g.Le Quéré et al., 2007;Lovenduski and Ito, 2009).Such a response on ocean carbon fluxes has been highlighted by Lovenduski et al. (2008) on shorter time scales.In this study, authors demonstrated that the southward shift of westerlies associated with the positive trend in the Southern Annular Mode (SAM) could drive a weakening of the Southern Ocean carbon sink because of changes in water mass overturning.However, it seems that such a response of both water mass overturning and ocean carbon sink may be overstated at the decadal time scale since no changes in Southern Ocean circulation have been detected from observations, while these data detect coherent warming and freshening trends in subsurface waters (Boening et al., 2008).Therefore, from these previous studies, it remains unclear if detected changes in ocean carbon sink within the last decades can be attributed to anthropogenic climate change or to natural climate variability that is a combination of an internal variability and a naturally forced variability (related to volcanoes and solar forcings).One large uncertainty arises from the lack of knowledge on ocean carbon flux variability at decadal time scales (Bates, 2007;Takahashi, 2009).
So far, these time scales have been an active fields of research in the ocean dynamics community (e.g.Latif et al., 2006, and references therein).Studies generally agree to the fact that sea-surface temperature and other ocean variables exhibit low-frequency modes of variability within mid-and high latitude oceans (e.g.Latif et al., 2006;Zwiers, 1996;Boer, 2000Boer, , 2004)).It is, thus, likely that the modes of variability of the ocean carbon fluxes could be inherited from low-frequency modes of variability like the Atlantic Multidecadal Oscillation (AMO; e.g.Bates, 2012;McKinley and Follows, 2004;McKinley et al., 2011), or the Pacific Decadal Oscillation (PDO; e.g., Valsala et al., 2012;Takahashi et al., 2006) in a similar way as their interannual variability is mainly associated with the tropical Pacific ocean variability caused by El Niño and La Niña (e.g.Chavez et al., 1999;Wang and Moore, 2012;McKinley and Follows, 2004).
However, studying low-frequency variations in ocean carbon fluxes within these oceanic regions is not straightforward.The major limitation comes from the temporal coverage of the observations or model reanalysis.This limitation can be illustrated by considering the North Atlantic for investigating multi-year variations in ocean carbon fluxes (Ullman et al., 2009;Löptien and Eden, 2010;Thomas et al., 2008).Here, large datasets (e.g.Gruber et al., 2002;Bakker et al., 2012) are available and ocean reanalysis are consistent over the 30 to 50 last years.Despite this large amount of information, an optimal spatial and temporal coverage of 30 to 50 yr provided by model reanalyses only allows studying 3 to 5 independent decades.Such limitation is, thus, even stronger for in situ observations.Observations are either combined and interpolated into large ocean regions in order to ensure a significant time sampling (McKinley et al., 2011), or merely considered in the terms of individual longterm stations (e.g.Gruber et al., 2002;Bates, 2007;Bates et al., 2012) that may be poorly representative of the basinscale dynamics (McKinley and Follows, 2004).
In this context, long model simulations (> 500 yr) are the only way to circumvent spatiotemporal sampling issues and have been demonstrated to include the minimum years of data to assess significantly variability at decadal time-scales (e.g.Boer, 2004).The usefulness of long model simulations for studying temporal variations of land and ocean carbon fluxes was first illustrated in Doney et al. (2006).Based on a 1000 yr fully coupled climate-carbon cycle simulation, the authors provide a good overview of the natural variability of the global carbon cycle in CCSM-1.4.However, in this study, both decadal and oceanic aspects receive less attention than the interannual variability of the terrestrial carbon cycle.
Here, as a first step to assess the internal decadal variability of the ocean carbon fluxes, we used an extended 1000 yr long preindustrial simulation performed with the Earth System Model IPSL-CM5A-LR.Using output of this model simulation, we performed several statistical time-series analyses to (i) track the low-frequency modes variability, (ii) locate the oceanic regions that contribute the most to these modes and (iii) identify the main drivers of the ocean carbon fluxes variability at decadal time scales.

Model description and pre-industrial simulation
This study exploits model output of IPSL-CM5A-LR (Dufresne et al., 2013), which contributes to the ongoing Coupled Model Intercomparison Project (CMIP5).IPSL-CM5A-LR combines the major components of the climate system: the atmosphere and land models are the atmospheric general circulation model LMDZ5A (Hourdin et al., 2013) and the land-surface model ORCHIDEE (Krinner et al., 2005).The atmospheric and land components use the same regular horizontal grid with 96 × 95 points, representing a resolution of 3.75 • × 1.875 • , while the atmosphere has 39 levels on the vertical.The oceanic component  (Madec, 2008).It offers a horizontal resolution of 2 • refined to 0.5 • in the tropics and 31 vertical levels.NEMOv3.2includes the sea ice model LIM2 (Fichefet and Maqueda, 1997), and the marine biogeochemistry model PISCES (Aumont and Bopp, 2006).PISCES simulates the biogeochemical cycles of carbon, oxygen and nutrients using 24 state variables.Macronutrients (nitrate, ammonium, phosphate and silicate) and the micronutrient iron limit phytoplankton growth and thus ensure a good representation of the high-nutrient low-chlorophyll regions (Aumont et al., 2003;Aumont and Bopp, 2006).Inorganic carbon pools are dissolved inorganic carbon, alkalinity and calcite.Total alkalinity includes contributions from carbonate, bicarbonate, borate, hydrogen and hydroxide ions (practical alkalinity).
For dissolved CO 2 and O 2 , air-sea exchange follows the quadratic wind-speed formulation (Wanninkhof, 1992).The core of this study is an extended preindustrial simulation of 1000 yr, in which the initial state of marine carbon cycle comes from a 3000 yr offline spinup plus 300 yr of online adjustment, while the dynamical components of IPSL-CM5A-LR have been spun-up for 600 yr (Dufresne et al., 2013).Spin-up strategy and model-data skill-assessment of (Dufresne et al., 2013).Spin-up strategy and model-data skill-assessment of IPSL-CM5A-LR's modern state of marine biogeochemistry are presented and discussed in Séférian et al. (2013).

Analytical methodology
The goals and methodology of this study are twofold.First, we aim at tracking and quantifying low-frequency oscillations of ocean carbon fluxes (here decadal to multi-decadal internal variability).For that purpose, we employ several statistical time-series analyses like the continuous wavelet transform (Torrence, 1998), the spectral density or the principal component analysis (PCA;Von Storch and Zwiers, 2002).
Secondly, we decipher which physical or biogeochemical drivers (here alkalinity, dissolved inorganic carbon, salinity and temperature) control the various modes of variability of ocean carbon fluxes.We employ a decomposition of sea-water pCO 2 (driving the ocean carbon fluxes) based on Takahashi et al. (1993) assuming that pCO 2 is a function of Alkalinity (Alk), dissolved inorganic carbon (DIC), salinity (S) and temperature (T ): pCO 2 = f g (Alk, DIC, T , S).The variance of a function of several random variables is given by expanding the function f (X 1 , X 2 , X 3 , . . ., X n ) (here pCO 2 ) in Taylor series around the mean values of X (here considered as the 1000 yr yearly-averaged climatology of Alk, DIC, S and T ): where σ 2 X i is the variance of the random variable X i and Cov(X i , X j ) is the covariance between the random variables X i and X j .
Based on this variance decomposition, we have conducted several offline computations of monthly-averaged sea-water pCO 2 and ocean carbon fluxes in which only one ocean carbon flux driver is allowed to vary in time while the others are fixed to their 1000 yr yearly-averaged climatology (Table 1).Online and offline computations of monthly-averaged seawater pCO 2 and ocean carbon fluxes (f gCO 2 ) show small differences in global average long-term mean and compare well in terms of variability (with a correlation of R ∼ 0.96).
Consequently, we assume that our variance decomposition strategy allows us to single out the variance contribution of Alk for pCO 2 -Alk or f gCO 2 -Alk, DIC for pCO 2 -DIC or f gCO 2 -DIC and so on.A similar approach was used in previous studies (e.g.Ullman et al., 2009;McKinley and Follows, 2004).
As described in Boer (2000Boer ( , 2004)), drifts in the variables matter for assessing low-frequency modes of variability.We have, therefore, removed drifts of ocean carbon fluxes and carbon-related fields, which have been estimated from linear least-square regression in function of time for each variable.Yet, in our case drifts in sea-water pCO   average.Regionally, the simulated oceanic sinks and sources of carbon agree with those estimated by inverse modelling with a global root-mean-square error of about 0.06 Pg C yr −1 (Fig. 1b).Only in the high latitude North West Pacific and the low latitude North Atlantic, simulated ocean carbon fluxes differ in sign from those estimated by inverse modelling.In the southern sub-polar regions, the model biases can be explained by a shift of outgasing structures probably associated with the equator shift of main wind stress structures (Marti et al., 2009): compared to the inverse modelling-based estimates, our model simulates a stronger outgasing of carbon in the southern sub-polar Atlantic than in the polar Southern Ocean.
If we consider the time-variability over 1000 yr of the globally integrated ocean carbon flux in the preindustrial simulation, we can track low-frequency oscillations from its year-to-year variability (Fig. 2a).The 10 and 20 yr timefiltered (i.e., running mean) variances amount, respectively, to 0.08 and 0.05 Pg C yr −1 representing 56 and 37 % of the interannual variance (0.14 Pg C yr −1 ).
The imprint of low-frequency modes of variability on the global carbon flux can be diagnosed from the autocorrelation function (ACF) of the global ocean carbon flux (Fig. 2b).The shape of the ACF does not correspond to a first order autoregressive process (AR1), but rather to a third order autoregressive process (AR3) according to the Yule-Walker estimation (Von Storch and Zwiers, 2002).The fact that the ACF of the global ocean carbon flux cannot be mimicked by an AR1 indicates that long-term memory processes likely drive ocean carbon fluxes.This implies that the ocean carbon flux behaves as other oceanic variables, like sea-surface temperature or sea-surface salinity (Frankignoul, 1985), but its low-frequency variability has a strong imprint over several years.
Low-frequency modes of the global ocean carbon flux variability range within a time-window of about 128 to 200 yr (centennial mode) and 18 to 50 yr (decadal mode) as shown on the wavelet time-frequency diagram (Fig. 1c and d).In comparison to the centennial mode, the decadal mode is not continuously detected.Nevertheless, the wavelet power R. Séférian et al.: Decadal variability of ocean CO 2 fluxes spectrum within the decadal time-window is significantly stronger than a red noise wavelet spectrum (at 95 % in red dots, at 90 % in blue dots, Fig. 1d).The fact that the decadal mode appears sporadically over the time scales underlines the complexity of studying a globally integrated quantity, because such a behaviour can be explained by a combination of several regional modes of variability either in phase or out of phase.As a consequence, in the next section, analyses will be conducted regionally.

Tracking the decadal mode of variability of regional ocean carbon fluxes
Modes of variability of one variable are an expression of its variance at several time scales.Here, the way we apprehend these different modes of variability, in a first approach, is to look at the variance of the time-filtered carbon fluxes in several oceanic regions, which correspond to those defined and used in Mikaloff Fletcher et al. (2007).The time filtering consists in a running mean of the regional carbon fluxes at 5, 10 and 20 yr (Fig. 3a).Except in the Arctic, the 1, 5 and 10 yr variances of the high latitude oceans are stronger compared to those found in mid-and low latitude regions.Yet, it is striking to see how the Southern Ocean variances (i.e., for all time filtering) dominate those of the other oceanic regions with a standard deviation of about 0.1 Pg C yr −1 .Such differences between the Southern Ocean and other oceanic regions are evidently fueled by differences in surface areas, which are much larger in the Southern Hemisphere regions than in the northern ones.As explained in Mikaloff Fletcher et al. (2007), region boundaries have been estimated using in situ measurements of pCO 2 (Takahashi et al., 2002) and a general ocean circulation model, and have been designed consistently with TRANSCOM regions (Gurney et al., 2002).That is, each region differs greatly in surface area compared to each other.Therefore, surface area weighing allows comparing the different regions independently of their design (Fig. 3b).Yet, variability of surface area weighted carbon fluxes is still much stronger in mid-and high latitude oceans than in the tropical ones.
In order to assess the relative importance of low-frequency modes of variability comparatively to interannual variability, we have employed a signal-to-noise ratio, the diagnosed potential predictability (DPP; Boer, 2004;Boer and Lambert, 2008;Persechino et al., 2013): where σ 2 is the ocean carbon flux 1 yr variance and σ 2 M (M = 5, 10 and 20 yr) is the variance of the time-filtered ocean carbon flux.The DPP metric (Fig. 3c) indicates that the 5 yr and the 10 yr time-filtered variances of ocean carbon fluxes within the North Atlantic, the North Pacific and the Southern Ocean account, respectively, for more than 30 % of the 1 yr variance.The Southern Ocean is the only oceanic region where multidecadal variability of ocean carbon fluxes (here represented by the 20 yr time-filtered variance) reaches up to 40 % of the 1 yr variance.These results are consistent with previous studies showing that low-frequency variability of various climate variables is mainly located in the high latitude oceans (Boer, 2000(Boer, , 2004;;Persechino et al., 2013;Mauget et al., 2011).
In particular, the three oceanic regions are characterised by well understood modes of variability in ocean dynamical fields: the AMO for the North Atlantic (Kerr, 2000;Enfield et al., 2001), the PDO for the North Pacific (Mantua and Hare, 2002) and the signature of the Southern Annular Mode on the Southern Ocean (Thompson and Wallace, 2000;Thompson et al., 2000).These various modes of variability have been suggested to have an influence on the variability of ocean carbon fields (Bates, 2012;McKinley and Follows, 2004;McKinley et al., 2011;Valsala et al., 2012;Takahashi et al., 2006;Lenton and Matear, 2007;Lovenduski et al., 2008).However, the fact that the signal-to-noise ratio of ocean carbon fluxes are comparable with those of other climate fields at regional scale with the same model (e.g., seasurface temperature, sea-surface salinity; Persechino et al., 2013) does not mean that the modes of variability or even the drivers are the same.
To address this issue, we have computed empirical orthogonal functions (EOFs) of ocean carbon fluxes for the North Atlantic (10 • -80 • N), The North Pacific (10 • N-70 • N) and the Southern Ocean (< 50 • S).We have selected the five leading principal components (PCs), which explain together a minimum of 80 % of the regional carbon flux total variance.The ratio between regional and global interannual variances explained by the leading EOFs amounts to 12, 20 and 26 % in the North Atlantic, the North Pacific and the Southern Ocean, respectively.As shown on Fig. 4, the spectral density of the five leading PCs of the regional carbon fluxes contains strong signals at decadal to multi-decadal time scales and displays large differences compared to a theoretical red-noise spectrum (AR1).The energy contained in the first mode of variability is always significantly stronger than the secondary modes except for the Southern Ocean for which the second mode of variability is stronger than the first mode at 5-10 yr.
Figure 4a shows a strong 20 yr mode of variability in the North Atlantic carbon flux, which is also found in several PCs (i.e., 1st, 3rd and 5th).The leading PC that accounts for 54 % of the total variance is the most energetic mode of variability.Compared to the leading PC, the second and the third PCs only account for 12 and 9 % of the total variance.The leading mode of the North Atlantic carbon flux variability has an amplitude of about 0.008 Pg C yr −1 that has been estimated from the standard deviation of the leading CP × EOF.Escudier et al. (2013)    SSS, SST and sea-ice cover with implications on mixed-layer depth anomalies in the sub-polar gyre and Nordic Seas.
In comparison to the North Atlantic, the various modes of variability found in the North Pacific and the Southern Ocean carbon fluxes differ at the decadal time scales.In the North Pacific, modes of variability range between 10 to 30 yr (Fig. 4b).The three leading PCs, which account for 76 % of the total variance, have a large part of their energy in a time-window of 10 to 30 yr.In this time-window, the leading PC explains 44 % of the total variance and presents modes of variability at 10 and 30 yr.The standard deviation of the leading CP × EOF associated with the multi-decadal modes of variability contribute to ∼ 0.03 Pg C yr −1 to the North Pacific carbon flux variability.
In the Southern Ocean, the leading PC (34 % of the total variance) sheds light on a prominent mode of variability at 50 yr, while a decadal mode of variability is found in the second PC that amounts to 28 % of the total variance (Fig. 4c).Among these different modes, variability is stronger for the first 50 yr mode than for the others.This mode has a magnitude of about 0.05 Pg C yr −1 .Figure 4c shows also a centennial mode of variability that appears in all of the PCs.However, this mode of variability cannot be thoroughly studied with a 1000 yr time-series.Similar results are found for regional-averaged surface pCO 2 (not shown).

Identifying drivers of low-frequency variability of ocean carbon fluxes
A further step is to identify drivers and examine their respective role in explaining the low-frequency variability of ocean carbon fluxes described in the previous sections.To find out what drivers contribute the most to the low-frequency variability of the ocean carbon fluxes, we use the variance decomposition detailed in Eq. ( 1) and Table 1.With this variance decomposition, we have conducted PCA analysis to identify correspondence in patterns and time variability between the fully-driven ocean carbon fluxes (fgCO 2 ) and fluxes driven only by Alk, DIC, S and T , named respectively f gCO 2 -Alk, fgCO 2 -DIC, fgCO 2 -SSS and f gCO 2 -SST.Figures 5 to 7 show the leading EOF for the North Atlantic, the North Pacific and the Southern Ocean carbon fluxes, while lagged correlations between the leading PC of ocean carbon fluxes and dynamical drivers, i.e., sea-surface temperature, mixed-layer depth (MLD) and sea level pressure (SLP) are presented in Fig. 8.These dynamical drivers have been chosen since they are proxies of climate modes of variability (i.e., the AMO, the PDO, the NAM and the SAM).AMO and PDO have been estimated from the leading EOF and PC of SST, while NAM and SAM are calculated from those of SLP.Their definition differ from their canonical definition, especially for the AMO and PDO indices, which are usually estimated from the 10 yr running mean of detrended SST north of the Equator over the Atlantic (Enfield et al., 2001;Kerr, 2000) and Pacific (Mantua and Hare, 2002), respectively.Here, our definition of AMO and PDO indices on the basis of PCA analysis show a significant correlation (> 0.42) with their canonical indice estimated from the same 1000 yr long preindustrial simulation.Yet, the computation of these indices is sensitive to the latitudinal boundaries.By computing these climate indices within 0 • -80 • N, comparison between our definition and their canonical definition gives higher correlations (> 0.6).That is, our definition of AMO and PDO mode of variability can be understood as a good approximation of their canonical definition.Same conclusions can be drawn with the NAM and SAM indices since correlations between our definitions and their canonical definition is high (R ∼ 0.7).
In the North Atlantic, f gCO 2 exhibits a horseshoe pattern, in addition to a dipole structure located close to the ocean convection sites of the Labrador and the Irminger Seas.Such features are also found in the leading EOFs of f gCO 2 -DIC and f gCO 2 -SST (Fig. 5c and e).The leading EOF of the f gCO 2 -SST exhibits an AMO-like pattern that is also highlighted in Ullman et al. (2009) using the MIT OGCM reanalysis over 1992-2006.In our model, the leading f gCO 2 -SST mode of variability displays a good match with the AMO-like pattern (Table 2).The respective correlations between PCs indicate a better match between time variability of f gCO 2 and that of f gCO 2 -SST (R ∼ 0.6) than that of f gCO 2 -DIC (R ∼ 0.3).Regarding f gCO 2 -Alk and f gCO 2 -SSS, the leading EOF spatial patterns and the respective PCs do not display correlations as strong as f gCO 2 -DIC and f gCO 2 -SST with the f gCO 2 (Fig. 5b and d).
Regarding variations between SST, MLD and SLP with those of f gCO 2 , lagged correlations indicate that the leading mode of f gCO 2 variability is strongly in phase with that of SST (Fig. 8).Significant correlations are also found between the leading PC of SLP and that of f gCO 2 , but can be understood as a feedback of the dynamical coupling occurring between the atmosphere and the ocean in the North Atlantic sector.Indeed, recent studies based on model simulations (Gastineau and Frankignoul, 2012) or observation-derived climate indices (Rossi et al., 2011;D'Aleo and Easterbrook, 2011;McCarthy et al., 2012) have argued that oceanic influence through SST variations lead the Northern Hemisphere climate (e.g., SLP) by a few years explaining the slightly lagged inverse relationship found between the NAM/NAO and the AMO climate indices.This seems to be the case in our model since cross-correlation between the leading PC of SLP and that of SST (not shown) displays strong correlation at lag −1/0 of about 0.64.Correlation between the leading PCs of f gCO 2 and MLD seems to indicate that f gCO 2 variations lead those of MLD.Such a result may be explained by the large regional boundaries that include several sites of deep convection that are either in phase or out of phase.Nonetheless, the presence of 10-yr-lagged oscillations in correlation coefficients between MLD and f gCO 2 (also in SST) can be explained by the 20 yr cycle occurring in the North Atlantic sector in this model (Escudier et al., 2013).
In the North Pacific, the leading EOF presents a strong positive pattern in the northern part of the domain (Fig. 6a), which is located at the same place as the strong North Pacific carbon outgasing (Fig. 1a).Such a pattern appears for f gCO 2 -DIC, but not for fgCO 2 -Alk, fgCO 2 -SSS and fgCO 2 -SST (Fig. 6).The leading EOF of the fgCO 2 -SST exhibits a PDO-like pattern and presents good agreement with the PDO diagnosed from the leading SST EOF (Table 2).This implies that PDO imprint a low-frequency signature on the ocean carbon fluxes (Valsala et al., 2012), which is perturbed by the contribution of the other drivers to the variability of the North Pacific carbon flux.The temporal correlations of the respective driver-related PCs with that of f gCO 2 demonstrate that variability driven by DIC mainly controls this mode of variability (R > 0.9).
Lagged correlations (Fig. 8) between leading PCs of f gCO 2 and dynamical drivers point toward a significant role of the mixed-layer depth (R ∼ 0.4 at lag 0) and the SLP (R ∼ 0.1 at lag 0).This supports the fact that both Ekmaninduced upwelling and winter deep-mixing events bring DIC-rich deep waters at surface contributing to the f gCO 2 variability.The strong lagged correlation at 2 yr found in SST and SLP leading PCs are related to the mid-latitudes wind regimes over the North Pacific and can be, at worst, considered as an artifact of the large regional boundaries we have chosen to define the North Pacific (10 • N-70 • N).They do not appear if one focuses above 40 • N, while the strong correlation between f gCO 2 and MLD and SLP are reinforced (R > 0.5 at lag 0).
The Southern Ocean differs from other oceanic regions due to its zonal structure and its ocean dynamical features like the Antarctic Circumpolar Current, the westerlies forcings and the near-shelf dense water mass formation.Such dynamical features influence the ocean carbon cycle and are mirrored on the leading EOF (Fig. 7a).A linear decomposition of drivers shows that leading EOFs of the different driver-related f gCO 2 present a good agreement with that of f gCO 2 (R ≥ 0.3).This strong spatial correlation is mainly associated with the strong spot of variability close to the Antarctic shelf.Yet, the temporal variability of the f gCO 2 -DIC strongly agrees with that of f gCO 2 (R > 0.9) compared to the variability of the other driver-related carbon fluxes (Fig. 7a).
To better understand Southern Ocean modes of variability, we have considered two sub-regions: one close to the Antarctic shelf and the sea-ice border (< 65 • S) and the second that extends from 65 • S to 45 • S, strongly influenced by the westerlies.In the wind-driven region, the leading EOF of f gCO 2 -DIC strongly matches with the fgCO 2 one (R ∼ 0.8).This result agrees with several studies, which have demonstrated that westerlies influence the vertical supply of DIC through the outcrop of deep water masses (Lenton and Matear, 2007;Lovenduski et al., 2008).In this sub-region, the prominent mode of variability is at 10 yr, and appears consistently within the second mode of variability of Southern Ocean carbon fluxes (Fig. 4c).
In the sea-ice driven region, PCA analysis shows instead that f gCO 2 appears to be driven rather by SSS than by DIC with a spatial correlation of about 0.7 and 0.5, respectively.Yet, temporal correlations indicate a better correspondence with DIC-driven variability (R ∼ 0.87) than the SSS-driven variability (R ∼ 0.2).
Similarly to the North Pacific, it seems that the MLD and the SLP variations drive the variability of the f gCO 2 in the Southern Ocean, not the SST.Two different processes take place: on the one hand, close to the sea-ice shelf, the deepconvection events occurring during winter mix subsurface DIC-rich waters with surface waters.This process is linked to sea-ice variability that also drives a fraction of the SSS spatiotemporal variability.On the other hand, in the open sea, Table 2. Temporal correlation (R 2 ) at lag 0 between the Atlantic Multidecadal Oscillation (AMO), the Pacific Decadal Oscillation (PDO), the Northern Annular Mode (NAM) and the Southern Annular Mode (SAM) with the regionally integrated carbon fluxes related to each drivers (i.e., Alk-f gCO 2 , DIC-f gCO 2 , SSS-fgCO 2 , SST-f gCO 2 ).Correlations between dynamical indices and regionally integrated carbon fluxes are done in the same region and are assessed with t test at 95 % of significance are mentioned in bold (following Bretherton et al., 1999).AMO and PDO indices have been estimated from the leading PC of SST over the North Atlantic (10 • N-80 • N) and the North Pacific (10 • N-70 • N), respectively.NAM and SAM indices have been computed from the leading PC of SLP over the North Atlantic (10 the wind forcing influences the subduction of intermediate and mode waters around Antarctica (Sallée et al., 2012) and Ekman-induced upwelling of deep-core waters.At decadal time scales, it is likely that the SAM oscillations could influence the Deacon cell through the Ekman-induced transport bringing in turn DIC-rich deep waters (Table 2).This may explain why we find a lag at 5-8 yr for the cross-correlation between SLP and f gCO 2 leading PCs.In the whole Southern Ocean region (< 50 • S), this mode linking SAM and ocean carbon fluxes, appears only as a second mode with a correlation of 0.7 at lag 0 indicating that this mechanism is of second order compared to the deep-mixing events linked to ocean-sea ice interactions.
In order to show how variability in ocean dynamics can interact with water mass properties in terms of dissolved carbon, we show the long-term mean concentrations and associated standard deviations of DIC, Alk and Alk-DIC at the maximum (winter) of the mixed-layer depth (Fig. 9).Alk-DIC serves as a good approximation for the carbonate ion concentration (Sarmiento and Gruber, 2006).This is a proxy for the buffer capacity of seawater and, hence, the chemical capacity for the ocean to take up atmospheric CO 2 .In regions of water mass formation (i.e., North Atlantic), the water mass content in DIC is set by air-sea interaction and biology and is, hence, sensitive to both variations in atmospheric forcings and variations in surface ocean properties (i.e., SST).In regions of water mass outcrops and transformations (i.e., North Pacific and Southern Ocean), water mass content in DIC is instead enriched by remineralisation processes occurring below the surface.Consequently, response of ocean carbon fluxes driven at first order by the difference in partial pressure of CO 2 between the atmosphere and the surface waters (which have been replenished in DIC due to DIC-rich deep water mass upwelling) overcomes the variability inherited from ocean thermodynamical variables like the SST and SSS.Such processes are highlighted in Fig. 10 showing, respectively, long-term mean state and variance profile of pCO 2 , Alk-pCO 2 , DIC-pCO 2 , SSS-pCO 2 and SST-pCO 2 .Long-term mean vertical profiles of pCO 2 exhibits large differences compared to pCO 2 -Alk, pCO 2 -DIC, pCO 2 -SSS and pCO 2 -SST at 500 m between the North Atlantic, the North Pacific and the Southern Ocean (Fig. 10ac).The pCO 2 long-term mean profile in the North Atlantic displays an important vertical gradient in the first 100 m and a well-mixed profile below (Fig. 10a).This is not the case for the North Pacific and the Southern Ocean pCO 2 profiles, which exhibit strong gradients from surface to 500 m (Fig. 10b and c).Within these two regions, the variability of pCO 2 at 500 m (here the standard deviation) is 5-to 6-fold stronger than the North Atlantic one.
Figure 10j-l present results of a multiple linear regression on the basis of Eq. ( 1).The statistical t test on a parameter slope (here Alk, DIC, S and T ) can be considered as an attribution test.That is, there is attribution of the variability driven by a given parameter to the pCO 2 variability if the slope is close to 1.There is no attribution if the slope is negative or if its confidence interval at 95 % (given by vertical lines on Fig. 10j-l) on of the slope includes 0. Consequently, our attribution test demonstrates that the variations (Fig. 10d-f) of pCO 2 at depth are mostly explained by variations in temperature in the North Atlantic with an attribution factor close to 1, while North Pacific and Southern Ocean pCO 2 variations at depth are set by variations in DIC with an attribution factor close to 1. Others variables like S and Alk have a poor contribution to the pCO 2 variability.

Conclusions
In this study, using a 1000 yr long preindustrial simulation performed with the Earth System Model IPSL-CM5A-LR, we have examined the decadal internal variability of ocean carbon fluxes as a part of their natural variability.Our results suggest that the low-frequency variability of ocean carbon fluxes (and sea-water pCO 2 ) emerges from year-to-year variability in all high latitude oceans.Low-frequency modes of variability display substantial differences in their periods of oscillation between the various oceanic basins.The North Atlantic carbon flux exhibits a strong 20 yr mode of oscillation while the North Pacific and the Southern Ocean carbon fluxes present several modes extending over a decade or more.The prominent modes of variability are a 10 yr and a 30 yr mode in the North Pacific and a 50 yr mode in the Southern Ocean.These modes of variability of carbon fluxes the three oceanic regions is set in part by the large-scale circulation of water masses and their biogeochemical properties.That is, in regions of dense water mass outcrops and transformations, variations in vertical supply of dissolved inorganic carbon (owing to Ekman-induced upwelling or deepmixing events) are of larger amplitude than the variations inherited from thermodynamical properties of surface waters (owing to ocean-atmosphere or ocean-sea ice interactions).This is not the case in regions of dense water mass formation.The fact that low-frequency variations in ocean carbon fluxes are simulated in our model within high latitude oceans is consistent with previous studies conducted on dynamical fields like sea surface temperature, surface air temperature or precipitation (e.g.Boer, 2004Boer, , 2000)).Of particular interest, a similar study with the same model (i.e., IPSL-CM5A-LR) by Persechino et al. (2013) demonstrates that decadal variations of sea surface temperature amount for 20 % (with a maximum of 50 %) of their interannual variations within the North Atlantic sector.Comparatively, ocean carbon fluxes exhibit low-frequency variations much stronger (∼ 20-40 %) than those of sea surface temperature.Such differences between low-frequency variations of ocean carbon fluxes and other dynamical fields are even stronger in the Southern Ocean, where Boer (2004) shows that multi-model zonal average variations of surface air temperature at decadal time scales only accounts for 10 % of the interannual variability.In our case, decadal variations of ocean carbon fluxes within Subpolar and Polar Regions of the Antarctic sector amounts to up to 25 %.Such features can be explained by the statistical properties of ocean carbon fluxes, which cannot be approximated with a first-order autoregressive process (like other dynamical variables mentioned above) indicating that longterm memory processes likely drive ocean carbon fluxes (and potentially other biogeochemical fluxes like the carbon export, the primary productivity and the remineralisation).
Also related to long-memory processes, all of these oceanic regions exhibit multi-centennial modes of variability (> 150 yr) that cannot be assessed with a 1000 yr preindustrial simulation.Similar multi-centennial modes of variability have also been revealed in sea-surface temperature, seasurface salinity, sea-ice volume and Atlantic meridional overturning for the GFDL-CM2.1 climate model (Delworth and Zeng, 2012).Interestingly, this study has hypothesised that these modes are controlled by interactions between ocean and sea-ice dynamics.It is, thus, likely that such interactions could impact air-sea carbon fluxes.
At decadal time scales, processes related to these lowfrequency modes of variability agree with process-based studies (Ullman et al., 2009;Corbière et al., 2007;Löptien and Eden, 2010;Thomas et al., 2008;Metzl et al., 2010;Lovenduski et al., 2008;Takahashi et al., 2006), which suggest that variations in air-sea fluxes (or sea-water partial pressure) of carbon are driven either by the sea surface temperature or by the vertical supply of dissolved inorganic carbon.However, our results seem to indicate that vertical supply of dissolved organic carbon related to the atmospheric modes of variability are of major importance in the North Pacific and the Southern Ocean, and not directly in the North Atlantic.Indeed, in our model, atmospheric modes of variability like the North Atlantic Oscillation weakly impact the ocean carbon fluxes (or very locally) as mentioned in Keller et al. (2012), while they can be projected on ocean-atmosphere coupled modes like the Atlantic Multi-decadal Oscillation.Nevertheless, we shall keep in mind that these findings are subject to criticism because they are related to the fact that we use only one Earth System Model.Consequently, we can wonder: (1) how do modes of variability of different Earth System Models compare to each other?(2) Is there any consistency between their drivers?Therefore, studies like the one conducted by Keller et al. (2012) are of great interest in order to understand how a given Earth System Model behaves in a multi-model perspective.
Yet, even in a multi-model perspective such as CMIP5, one of the major uncertainties in the quantification of the internal variability of ocean carbon fluxes (as others variables) would be related to the horizontal resolution of the considered ocean model (> 1 • in most of the CMIP5 models).Recent literature (e.g.Penduff et al., 2010Penduff et al., , 2011;;Lovenduski et al., 2013) highlights that meso-scale activity in eddying-regions like the Southern Ocean strongly enhances interannual variability of ocean variables -such as ocean turbulence.The impact of meso-scale activity could likely amplify the interannual variability of the ocean carbon fluxes.Since the variability of the ocean carbon fluxes is influenced by the long-term mean state of ocean carbon-related fields, this variability could be affected by the response of westerlies-induced transport to the model resolution (e.g.Farneti et al., 2010;Hofmann and Morales Maqueda, 2011) and, hence, model biases (Swart and Fyfe, 2011).
Notwithstanding the limitations of the model, our findings provide a first step to better understand the role of the internal variability (as a part of the natural variability) versus that of the anthropogenic forced variability on the ocean carbon fluxes.Our results demonstrate that care should be taken while analysing short-term changes with delta or biases correction methods, which consist in applying model anomaly to observed fields (e.g.Sarmiento et al., 1992).Indeed, these methods generally assume that the long-term mean state does not affect the variability, while our results demonstrate that it does matter in some oceanic regions like the Southern Ocean and the North Pacific.This is not the case for long-term changes for which rising atmospheric CO 2 and climate change can impact ocean carbon fluxes much more strongly than their decadal variations (Matear and Lenton, 2008;Lenton and Matear, 2007;Séférian et al., 2012;Roy et al., 2011).
The issue of what observations may be most useful to constrain predictions of the ocean carbon fluxes within decades has received much less attention than the issue of disentangling the decadal natural variability from its secular anthropogenic trend (McKinley et al., 2011).Both issues are important, but the observational needs are different.It seems that the observational needs are larger in oceanic regions where low frequency modes of variability take place than in those dominated by interannual variability.Yet, relevant spatiotemporal scales are unknown to ensure an optimal and efficient sampling strategy.To assess and quantify the observational need in the view of studying the natural variability of ocean carbon fluxes, similar studies to this of Lenton et al. (2009) could be a possible alternative.
2 and ocean carbon fluxes are small compared to the interannual or lowfrequency oscillations: for global ocean carbon fluxes it represents about −0.03 Pg C over 100 yr in the first century of the 1000 yr simulation and only −0.01 Pg C over 100 yr in the last century of the 1000 yr simulation.
have found a similar signature of a 20 yr mode in the North Atlantic in the same model and have explained this mode of variability by the interplay between

Fig. 5 .
Fig. 5. Spatial pattern of the leading empirical orthogonal function (EOF, in normalised unit) of (a) the fully-driven carbon fluxes, and that of the Alk-driven (b), DIC-driven (c), SSS-driven (d) and SST-driven (e) ocean carbon fluxes in the North Atlantic.Spatial correlations (based on EOFs) and temporal correlations (based on CPs) between the fully-driven ocean carbon fluxes and each driver-related carbon fluxes are mentioned in brackets.

Fig. 6 .
Fig. 6.Spatial pattern of the leading empirical orthogonal function (EOF, in normalised unit) of (a) the fully-driven carbon fluxes, and that of the Alk-driven (b), DIC-driven (c), SSS-driven (d) and SST-driven (e) ocean carbon fluxes in the North Pacific.Spatial correlations (based on EOFs) and temporal correlations (based on CPs) between each driver-related and the fully-driven ocean carbon fluxes are mentioned in brackets.

Fig. 7 .
Fig. 7. Spatial pattern of the leading empirical orthogonal function (EOF, in normalised unit) of (a) the fully-driven carbon fluxes, and that of the Alk-driven (b), DIC-driven (c), SSS-driven (d) and SST-driven (e) ocean carbon fluxes in the Southern ocean.Spatial correlations (based on EOFs) and temporal correlations (based on CPs) between the fully-driven ocean carbon fluxes and each driver-related carbon fluxes are mentioned in brackets.

Fig. 8 .
Fig.8.Lagged correlation of leading principal components of ocean carbon fluxes and sea-surface temperature (SST), mixed-layer depth (MLD) and sea-level pressure (SLP) in the North Atlantic, the North Pacific and the Southern Ocean.Null hypothesis assessed with a t test at 95 % of significance are mentioned with dashed blue lines and red lines (followingBretherton et al., 1999).

Fig. 10 .
Fig. 10.Scaled spectral density of the five leading principal components of ocean carbon fluxes for (a) the North Atlantic (10 • N-80 • N), (b) the North Pacific (10 • N-70 • N) and (c) the Southern ocean(< 50 • S).Dashed lines represent fitted AR1 theoretical spectrum according to the Yule-Walker estimation, while the linear dashed lines indicate the maximum of spectral density for each principal component.Confidence intervals at 95 % of significance of the scaled spectral density maximum are indicated by vertical bars; they have been computed from a (scaled) χ 2 distribution.

Table 1 .
Description of the distinct ocean carbon fluxes (fgCO 2 ) or oceanic partial pressure of carbon (pCO 2 ) driven by the variability of solely one driver (i.e., T/SST, S/SSS, DIC or Alk) compared to the fully-driven ocean carbon fluxes.