Variability of global net sea–air CO2 fluxes over the last three decades using empirical relationships

The interannual variability of net sea–air CO2 flux for the period 1982–2007 is obtained from a diagnostic model using empirical subannual relationships between climatological CO2 partial pressure in surface seawater (pCO2SW) and sea surface temperature (SST), along with interannual changes in SST and wind speed. These optimum subannual relationships show significantly better correlation between pCO2SW and SST than the previous relationships using fixed monthly boundaries.Our diagnostic model yields an interannual variability of±0.14 PgC yr−1 (1σ)with a 26-year mean of −1.48 PgC yr−1. The greatest interannual variability is found in the Equatorial Pacific, and significant variability is also found at northern and southern high-latitudes, depending in part, on which wind product is used. We provide an assessment of our approach by applying it to pCO2SW and SST output from a prognostic global biogeochemical ocean model. Our diagnostic approach applied to this model output shows reasonable agreement with the prognostic model net sea–air CO2 fluxes in terms of magnitude and phase of variability, suggesting that our diagnostic approach can capture much of the observed variability on regional to global scale. A notable exception is that our approach shows significantly less variability than the prognostic model in the Southern Ocean.


Introduction
The estimated annual CO 2 emissions attributable to anthropogenic activities are not balanced, on a year to year basis, with the rates of atmospheric CO 2 increase and long-term estimates of net CO 2 uptake of land and ocean (Sarmiento and Gruber, 2002). The annual imbalance estimated for each year since 1960 ranges from +3 to −2 PgC yr −1 (1 PgC = 10 15 g C) (Le Quéré et al., 2009). These imbalances are as large as the magnitude of annual ocean uptake rate of 2 PgC yr −1 , and may be attributed partially to the errors in the estimated magnitude and variability of the CO 2 sinks as well as to differences in response time of each CO 2 reservoir. Three-dimensional (3-D) global biogeochemical ocean circulation models (GCMs) generally yield small interannual variability ranging from ±0.2 to 0.35 PgC yr −1 (1σ ) (e.g. Le Quéré et al., 2003;Obata and Kitamura, 2003;McKinley et al., 2004;Doney et al., 2009a).
The partitioning of CO 2 uptake change between the land and the ocean has been determined with atmospheric records of CO 2 concentration, δ 13 CO 2 , and O 2 /N 2 (e.g. Francey et al., 1995;Keeling et al., 1995;Bender et al., 2005;Patra et al., 2006;Rayner et al., 2008). In these studies, the combined CO 2 uptake of land and ocean calculated from the difference between total CO 2 emissions and the observed increase in atmospheric CO 2 concentrations are separated by changes in atmospheric δ 13 CO 2 and O 2 /N 2 as proxies for CO 2 uptake by land biosphere. In contrast to the GCMs, these top-down approaches show relatively high magnitude of interannual changes in the ocean carbon reservoir. Alden et al. (this issue) suggest that the large variability of the ocean CO 2 sink obtained from atmospheric inversion method using δ 13 CO 2 could be caused by uncertainty in isotope disequilibrium fluxes and poorly quantified variability of isotopic fractionation in land plants. A joint atmosphere-ocean inversion method using atmospheric and oceanic observational constraints has been developed to more accurately estimate CO 2 fluxes in the Earth System . Results from recent 3-D atmospheric inversion models with long-term atmospheric sampling networks (Bousquet et al., 2000;Rödenbeck et al., 2003) suggest lower interannual variability in sea-air CO 2 exchange, more in agreement with those predicted by empirical methods (Lee et al., 1998;Park et al., 2006) and GCMs (e.g. Le Quéré et al., , 2003McKinley et al., 2004;Wetzel et al., 2005;Doney et al., 2009a). However, empirical approaches and ocean GCMs have shown interannual variability at the lower end of the reported range (±0.2 to 0.5 PgC yr −1 ). Key unresolved issues are which regions significantly contribute to the global ocean flux variability and how large-scale climate reorganizations such as the El Niño-Southern Oscillation (ENSO) and the North Atlantic Oscillation (NAO) affect variability (Le Quéré et al., 2003;McKinley et al., 2004McKinley et al., , 2006Peylin et al., 2005). The ocean carbon database is too sparse, except for a few regions, to quantify sea-air CO 2 flux variability on the required spatial and temporal scales from CO 2 observations alone, even with the recent increases of CO 2 partial pressure (pCO 2 ) measurements on moorings and ships of opportunity (SOOP). Lee et al. (1998) developed a diagnostic model based on the Takahashi et al. (1997) monthly climatology of surface water and atmospheric partial pressure differences ( pCO 2 = pCO 2SW − pCO 2AIR ) and empirical algorithms between pCO 2SW and sea surface temperature (SST), and inferred an interannual variability of ±0.2 PgC yr −1 (1σ ). Park et al. (2006) presented a more detailed description of this diagnostic model and evaluated uncertainties with time-series observations. This work is a follow-up to the aforementioned publications by Lee et al. (1998) and Park et al. (2006) with the following important improvements. The work is extended to cover 26 years for which consistent high-resolution SST and wind products are available. The new sea-air CO 2 flux climatology of Takahashi et al. (2007Takahashi et al. ( , 2009a is used as the basis for the optimum subannual pCO 2SW -SST relationships. The subannual pCO 2SW -SST relationships have been changed from three fixed time period relationships, spanning January-April, May-August and September-December, to grid cell specific variable time period relationships to define the subannual pCO 2SW and SST trends. The effect of wind speed products on CO 2 flux variability is also considered. There remain several assumptions in this method that cannot be validated on first principles. This approach assumes that subannual pCO 2SW -SST relationships derived from the pCO 2 climatology for each 4 • × 5 • grid cell can be applied to interannual SST anomalies to predict interannual variations in pCO 2SW . This means that for every grid cell variations in pCO 2SW can be captured with changes in SST and that these derived suban-nual linear trends are also applicable to determine pCO 2SW on interannual time scales. This diagnostic model cannot capture long-term trends unrelated to trends in SST such as those caused by increasing rate of release of fossil-fuel CO 2 . The assumptions are not easy to verify over the global ocean due to a dearth of measurements. The comparison presented here with two ocean time-series stations and particularly from applying the approach to the output of a GCM suggests that the approach captures a significant fraction (≈70%) of the variability in sea-air CO 2 fluxes.
A validation is performed by applying our approach to the output of a GCM. The model output of pCO 2SW is the result of the major biogeochemical processes and transport that occur in the ocean (Doney et al., 2009a). Optimum subannual relationships are created from the model pCO 2SW and SST in each grid cell in the same manner as was done with the pCO 2 climatology of Takahashi et al. (2009a). They are then applied to the SST anomalies produced by the model. This provides an important diagnostic on how well our method reproduces the interannual changes in pCO 2 on global scale. Much of the regional interannual variability determined with our approach appears related to the large-scale climate reorganizations associated with the ENSO in the Equatorial Pacific and the NAO in the North Atlantic.

Calculation method
The monthly mean net sea-air CO 2 flux (F ym ) for each 4 • (latitude) × 5 • (longitude) grid cell for an individual year other than 2000 is estimated from the global pCO 2 climatology and SST anomalies compared to the SST for reference year 2000 ( SST ym-2000m ) using subannual pCO 2SW -SST relationships [(δpCO 2SW /δSST) 2000m ] and gas transfer velocity in the following manner: or where subscript ym is the year (y) and month (m) during the time period of 1982-2007, and subscript 2000m refers to the month in 2000. The determination of (δpCO 2SW /δSST) 2000m is described in detail later. The solubility of CO 2 , K 0ym , is determined from monthly SST and monthly climatological sea surface salinity (SSS) estimates using the solubility equations of Weiss (1974 Takahashi et al. (2009a). Positive F ym values indicate that CO 2 is emitted from the ocean, whereas negative values indicate that the ocean is a CO 2 sink. The monthly gas transfer velocity, k ym for each grid cell is estimated from the second moment of monthly mean wind speed and the gas transfer coefficient where <U 10ym 2 > is the second moment of the wind at 10 m above sea surface representing the variance of the 6-hourly wind speeds for each grid cell over ym, and Sc is Schmidt number. The proportionality coefficient of 0.22 for <U 10 2 > is derived from the coefficient of 0.26 for monthly mean wind speed (<U 10 >) and global mean <U 10 2 >/<U 10 > 2 of 1.2 for ice-free oceans (0.26/1.2 = 0.22) used by Takahashi et al. (2009a). The coefficient of 0.26 differs from that of 0.39 proposed by Wanninkhof (1992) and is based on an updated global gas transfer velocity based on the partitioning of the global bomb-14 C inventory between atmosphere and ocean utilizing a global ocean circulation model (Sweeney et al., 2007). A different wind speed product is used as well in determining the coefficient compared to the original estimate of Wanninkhof (1992). For global gas exchange-wind speed relationships the gas transfer coefficient is intrinsically tied to the wind speed (Naegler et al., 2006;Sweeney et al., 2007), and 0.26 is appropriate for the winds used here. The monthly mean second moments are from the 6-hour National Centers for Environmental Prediction/Department of Energy (NCEP/DOE) reanalysis 2 product on a Gaussian grid (http://www.cdc.noaa.gov/ data/gridded/data.ncep.reanalysis2.html), henceforth called NCEP-II.
In the polar regions where sea-ice forms seasonally, we correct the sea-air CO 2 flux by assuming no flux through sea-ice. The monthly fractional sea-ice cover values for each 4 • × 5 • grid cell are obtained from the NCEP/DOE reanalysis 2 surface ice concentration fields (ftp://ftp.cdc.noaa.gov/Datasets/ ncep.reanalysis2/gaussian_grid/). The original data are regridded to a 4 • × 5 • grid and averaged for each month in each grid cell. Following the convention in Takahashi et al. (2009a), each grid cell is regarded as a sea-ice-free area when the ice cover value is less than 0.1. In the case that the ice cover value is over 0.9, we assume that each grid cell has 10% ice-free open water (ice cover value = 0.9) because of leads and polynyas where CO 2 is exchanged across the sea-air interface (Takahashi et al., 2009a).
As shown in eq. 2, the net sea-air CO 2 flux can be expressed in terms of pCO 2 and the variation of pCO 2SW arising from the SST anomaly relative to 2000. This means that the pCO 2 fields of Takahashi et al. (2009a,b) can be used for this work. The approach implicitly assumes that pCO 2SW increases at the same rate as the atmosphere over the time, modulated by SST anomalies. Takahashi et al. (2009a) show that over broad regions, this is a reasonable assumption but regional deviations in the increase in pCO 2SW not closely correlated with SST have been observed (see, e.g. Schuster and Watson, 2007;Metzl, 2009 fig. 1a in Takahashi et al. (2009a)]. The database used in this climatology is about three times larger than the previous climatology (Takahashi et al., 2002). The increased data allow for a better resolution of the seasonal trends and are believed to have fewer artefacts caused by interpolation. This allows for a more reliable estimate of the base state and determination of the subannual pCO 2SW -SST relationships. We use the version of the pCO 2 climatology provided in the website listed earlier. This is an update of the data used in Takahashi et al. (2009a) corrected for some minor errors.
Climatological SST values (monthly mean NOAA OI SST for the period of 1981-2006 without the El Niño months) used here for deriving subannual pCO 2SW -SST relationships differ from those based on the SST measurements concurrently made with the pCO 2SW in the climatology (Takahashi et al., 2009a). This is due to the relative sparsity of observational SST data in each grid cell of the pCO 2SW climatology compared to the NOAA OI SST product and use of different smoothing and interpolation functions in the NOAA OI SST product. Although these two sets of SST are closely in agreement for most of the cells, the SSTs in some cells, particularly in regions with high SST gradients differ as much as 7 • C. Nevertheless, the mean difference between the concurrent SST measurements and the NOAA OI SST product is small (+0.08 ± 1.4 • C, n = 21,108; Takahashi et al., 2009a), and hence the more complete NOAA OI SST products may be treated as consistent with the pCO 2 values in most of the grid cells.

Subannual relationships between pCO 2SW and SST, (δpCO 2SW /δSST) 2000m
A critical part of the analysis is the determination of the trends of pCO 2SW and SST for each grid cell of the monthly climatology. To avoid spurious values of (δpCO 2SW /δSST) 2000m that can arise by determining the trend from month to month, we determine a linear relationship of pCO 2SW and SST over seasonal or longer time periods, and use the slopes of these linear segmented relationships along with the SST anomaly ( SST ym-2000m ), for the months that are included in each relationship, to determine the anomaly in pCO 2SW (see eqs 1 and 2). As shown in the examples in Fig. 1, the relationships are discontinuous, leading to, sometimes, abrupt changes in seasonal response of pCO 2SW to SST anomalies. The relationships between pCO 2SW and SST for each grid cell are different than used in Lee et al. (1998) and Park et al. (2006). The new method called the optimum subannual pCO 2SW -SST approach is applied to each 4 • × 5 • grid cell, except for the central and eastern Equatorial Pacific. In this method we do not use fixed monthly boundaries as used in the previous analyses but rather pick boundaries based on an optimum segmented linear fit to pCO 2SW and SST data. Least squares linear fits of monthly pCO 2SW and SST values are determined for at least 3 consecutive months according to the annual patterns of pCO 2SW and SST (Figs 1e-h). The subannual relationships are created from 3 to 12 consecutive months, which are determined in the following iterative fashion (Park et al., 2010). First, for each grid cell, a linear fit is determined using 12 months of data. If the resulting correlation coefficient (R 2 ) is over 0.9, we accept this single pCO 2SW -SST relationship. If not, we determine subannual trends using time windows of at least 3 consecutive months. We use the following criteria to choose the time windows and create the corresponding set of optimum relationships for each grid cell: (1) all correlation coefficients exceed 0.5; (2) maximum mean value of correlation coefficients; (3) minimum number of relationships and maximum number of months in each relationship. When relationships have similar correlation coefficients (3) takes precedence over (2). If for a particular grid cell we cannot find a set of subannual relationships where all R 2 values exceed 0.5, we select subannual boundaries such that more relationships have R 2 greater than 0.5 and the set has a higher mean value of R 2 . Although the selections are somewhat subjective, the comparison below with fixed boundaries shows the significantly improved fits in the approach used here.  derived from multiyear observations that were collected on board SOOP and research ships from 1979 to 2008 (R. A. Feely, unpublished data). The drivers of interannual variability in this region are different from those of seasonal variability, and extensive pCO 2SW observations are available for the period of investigation covering seven El Niño and five La Niña periods. Therefore, unique pCO 2SW -SST equations for three different time periods (1979-1989, 1990-mid-1998 and mid-1998-2008) are used in this study. These equations are updated and extended through 2008 from those of Feely et al. (2006). Mean atmospheric pCO 2 values for estimating pCO 2 in each grid cell of the central and eastern Equatorial Pacific are calculated from mole fractions of atmospheric CO 2 obtained from the National Oceanic and Atmospheric Administration/Earth System Research Laboratory (NOAA/ESRL) (GLOBALVIEW-CO2, 2009).

Rationale behind the pCO 2SW -SST relationships
Surface pCO 2SW changes with physical chemical processes (changes in temperature, salinity and inorganic carbon speciation), biological process (photosynthesis and oxidation), transport process (mixing and advection of water masses with different CO 2 concentrations) and sea-air CO 2 flux (e.g. Takahashi et al., 2002). These processes are often related to changes in SST and this feature is used to create the empirical relationships we use here. Although the pCO 2SW shows strong seasonal and regional correlations with SST over much of the world's ocean, their patterns differ for different regions and seasons (e.g. Takahashi et al., 1993;Lefèvre and Taylor, 2002;Feely et al., 2006). For example, in subtropical oceans where biological production is low, the pCO 2SW is primarily regulated by SST, and hence increases as water warms from winter to summer, exhibiting a positive slope in the pCO 2SW -SST regression (Figs 1e and S1). In contrast, in subpolar oceans, deep waters rich in CO 2 mix upward to the surface regime as a result of winter cooling of surface waters, thus causing high pCO 2SW in spite of colder SST. In the spring, pCO 2SW is reduced due to photosynthetic utilization of CO 2 , although spring warming counteracts the biological effects partially. Thus, pCO 2SW tends to decrease with increasing SST, exhibiting a negative slope as seen in examples in Figs 1g and h.
In the previous studies, the seasonal pCO 2SW -SST relationships were determined from monthly pCO 2SW and SST values for three fixed time periods: January through April (Season 1), May through August (Season 2) and September through December (Season 3) (Lee et al., 1998;Park et al., 2006). Low correlations between pCO 2SW and SST were found in several regions and seasons in those studies particularly in the South Indian and Southern oceans. When three fixed seasonal relationships are derived from the updated monthly pCO 2SW climatology (Takahashi et al., 2009a), 62% of total grid cells for Season 1 show correlation coefficients less than 0.5 (Figs 1c and d). When the optimum subannual approach is used for each grid cell, only 1% of total grid cells show correlation coefficients less than 0.5 between SST and pCO 2SW (Fig. S2). The mean correlation coefficient for all the grid cells and all subannual periods is 0.83 ± 0.14 indicating that the low correlations in the previous approach are largely due to applying the criterion of fixed monthly boundaries. This suggests that SST change can reasonably characterize pCO 2SW variability in most grid cells of the climatology on subannual scales.
Examples of comparisons between three fixed seasonal relationships and the optimum subannual relationships for specific grid cells are illustrated in Fig. 1. For the optimum subannual relationships, about half of the total grid cells have three subannual relationships, but only 20% of them have the same monthly ranges as the previous seasonal relationships with the fixed monthly boundaries. Twenty-four percent of the grid cells have two subannual relationships. Eighteen percent of the grid cells have a single pCO 2SW -SST relationship, mainly in the northern subtropical regions. Only 7% of total grid cells have four subannual relationships. At high latitude, 10 grid cells show no meaningful pCO 2SW -SST relationships because there are no appreciable temporal variations of climatological pCO 2SW or SST. Global distribution of optimum subannual relationships for each month is shown in Supporting Information (Fig. S1).
The effect of this updated procedure for creating subannual pCO 2SW -SST relationships on pCO 2 estimation differs for each grid cell. The grid cells that have low correlation coefficients for relationships with the fixed monthly boundaries and with larger interannual SST variations show the greatest changes of interannual variability in annual mean pCO 2 (Fig. 2). Fiftyfive percent of total grid cells show increase of interannual variability in net sea-air CO 2 fluxes and 30% of them have lower variability compared to the old procedure using fixed relationships (Fig. S3). However, the interannual variability of global sea-air CO 2 fluxes estimated from the two relationships are the same (0.14 PgC yr −1 ; 1σ of annual global net sea-air CO 2 fluxes). This is because in the new scheme there is more interregional compensation that damps the global interannual variability. CO 2 fluxes estimated from optimum subannual relationships show 65% larger variability in the northern high-latitude and 15% larger variability in the subtropics compared to fluxes obtained from relationships with the fixed monthly boundaries (Table S1).

Effect of wind speed product
Wind speed directly controls the magnitude of sea-air CO 2 flux for a given pCO 2 (eq. 1), and temporal and spatial variabilities in wind speed contribute to variability in net sea-air CO 2 fluxes. To estimate the impact of different wind speed products on the interannual variability of CO 2 fluxes, we compare results using the NCEP-II product with results using a new cross-calibrated, multiplatform (CCMP) ocean surface wind 1982 1986 1990 1994 1998 2002 2006 ΔpCO 2 anomaly (μatm) Year 1982 1986 1990 1994 1998 2002 2006 ΔpCO 2 anomaly (μatm) velocity product for the 20-year period from 1988 to 2007 (Ardizzone et al., 2009). The CCMP wind speed product is based on multiple input data sources including three types of microwave radiometer sensors and scatterometers, to generate a wind speed product with less spatial smoothing (see http://podaac.jpl.nasa.gov/DATA_CATALOG/ccmpinfo.html).
The <U 10 > and <U 10 2 > values are 8.1 m s −1 and 83.8 m 2 s −2 for the NCEP-II wind and 7.4 m s −1 and 67.2 m 2 s −2 , for the CCMP wind, and the products show different regional distributions (Also see, Wallcraft et al., 2009 for a comparison of satellite winds and winds obtained from numerical weather products). Because the coefficient of proportionality in the gas transfer velocity-wind speed relationship 'a', k = a<U 10 2 >, obtained from global 14 C ocean inventory depends on wind speed product (Naegler et al., 2006), we recalculate a coefficient value of 0.27 for the second moment of CCMP wind speeds. This coefficient is obtained from the ratio of the global average second moments of the CCMP and NCEP-II products multiplied by the proportionality coefficient of 0.22 used with the NCEP-II wind product.
Global mean sea-air CO 2 flux values for the period over which both wind speed products are available (1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007) are −1.49 ± 0.14 PgC yr −1 (1σ ) for the NCEP-II wind product and −1.29 ± 0.14 PgC yr −1 (1σ ) for the CCMP wind product (Fig. 3a). Thus, even after a 20% reduction in the proportionality coefficient to account for differences in global mean second moment of wind speeds, the NCEP-II wind speed product yields a global net sea-air CO 2 flux 16% higher than the CCMP wind speed product. This is attributed to regional differences in the wind products. The CCMP wind speeds are lower than the NCEP-II winds at mid-and high-latitudes that are predominantly CO 2 sink regions. The CCMP winds are higher in the tropical areas, which are generally sources of CO 2 to the atmosphere. Higher CO 2 efflux in the tropics and lower oceanic CO 2 influx in the higher latitudes using the CCMP wind product lead to the lower global net sea-air CO 2 flux (Fig. 3b).
The variability relative to the global mean CO 2 flux using the CCMP wind product is 20% higher. The regional patterns of CO 2 flux variability are also different (Fig. 3c). Predominant differences are found in the Equatorial Pacific, and differences are also observed in the mid-latitude and northern high-latitude oceans. Selection of wind speed product significantly impacts not only the estimate of the global sea-air CO 2 flux but also its interannual variability.

Interannual variability of global net sea-air CO 2 fluxes
The 26-year mean regional CO 2 fluxes and their variability are presented in Table 1. The variability is expressed as a standard deviation (1σ ) of annual net sea-air CO 2 fluxes. The Equatorial Pacific shows the largest interannual CO 2 flux variability (18%) relative to the mean efflux 0.44 PgC yr −1 for the period 1982-2007. The Southern Ocean and regions at northern high-latitudes show significant variability of 11% and 10%, respectively, compared to the 26-year mean. The global net sea-air CO 2 fluxes show a strong correlation with the ENSO cycle. Higher oceanic CO 2 uptakes (negative anomalies) occur during the El Niño periods (Fig. 4). Our analysis shows that the interannual variability of global net sea-air CO 2 flux is ±0.14 PgC yr −1 (1σ ) for the period 1982-2007. This value is smaller than the variability previously estimated in Lee et al. (1998) and Park et al. (2006), largely due to the difference of proportionality coefficient used in gas transfer velocity calculation and the different wind speeds (U 10 ). The relative interannual variations compared to the global mean net sea-air CO 2 flux values are comparable since the net global uptake in the previous works is larger as well. High-North, northern high-latitude oceans poleward of 42 • N, and Southern Ocean (>42 • S). c The percentage of CO 2 variability relative to the regional 26-year mean annual CO 2 flux.  (Trenberth, 1997). For comparison with ONI, annual net sea-air CO 2 flux values are plotted in the middle of each corresponding year. The 26-year mean sea-air CO 2 flux value and variability for each region are presented in Table 1.
Tellus 62B (2010), 5 Regionally, the eastern Equatorial Pacific and the SE Pacific off the western coast of South America show large interannual variability closely connected to the ENSO cycle as detailed later (Figs 4 and 5). The high-latitude North Atlantic (north of 50 • N; Fig. 5), a region of large CO 2 uptake, also shows large interannual variability due to large SST anomalies and strong subannual trends between pCO 2SW and SST. The Indian Ocean shows relatively low interannual variability despite interannual changes in the seasonal upwelling due to year-to-year changes in Southwest Monsoon winds. This is because the region directly impacted by the strong Southwest Monsoon winds is relatively small, and there is little variability in non-upwelling seasons and regions. The sea-air CO 2 flux in the Southern Ocean is a function of significant drawdown of surface pCO 2SW by biological productivity in summer and release of CO 2 due to upwelling of deep water rich in CO 2 in winter. Grid cells showing high interannual variability near 60 • S (red colour in Fig. 5b) correspond to the grid cells located along the edge of the seasonal sea-ice field that experience CO 2 effluxes in August [see fig. 15b in Takahashi et al. (2009a)]. The pCO 2SW -SST relationships are strongly negative, which combined with large interannual SST anomalies leads to large interannual variability in sea-air CO 2 fluxes.
Long-term trends in regional sea-air CO 2 fluxes are suggested by our diagnostic approach (Fig. 4). These trends are caused by long-term trends in SST, because we assume that pCO 2 only changes by surface ocean processes related to SST variations (eq. 2). Increase in oceanic CO 2 uptake in northern high-latitude oceans of −0.05 PgC decade −1 (grey colour in Fig. 4) is related to long-term increase of SST in the North Atlantic (>42 • N) and predominant negative pCO 2SW -SST relationships in this region. Recent observational and model studies highlight that trends of CO 2 fluxes are caused by climate-induced physical changes not necessarily correlated with SST (Le Quéré et al., 2007;Schuster and Watson, 2007;Lovenduski et al., 2008). Thus, our approach for diagnosing interannual variability of CO 2 fluxes will not adequately capture decadal and much longer time scale trends of CO 2 fluxes. Even with long-term trends in SST, it is not clear that subannual regressions are appropriate estimates of the trend in pCO 2 .

Correlation between CO 2 fluxes and large-scale climate reorganizations
Changes in natural modes of atmosphere-ocean coupled variability, such as the ENSO, NAO and Southern Annular Mode (SAM) cause large variations in climate and weather over much of the globe on interannual and longer time scales. They also have striking impact on oceanic variability through associated changes in heat content, water circulation and winds that could affect pCO 2 and sea-air CO 2 exchange. The impact of these climate modes on sea-air CO 2 fluxes is of particular interest to determine if our approach of using subannual regressions of pCO 2SW and SST to estimate interannual variability can capture changes in sea-air CO 2 fluxes caused by these multi-annual climate modes.

The El Niño-Southern Oscillation.
ENSO is an oscillation of the ocean-atmosphere system in the tropical Pacific that affects weather around the globe. It is often characterized by changes of SST in the eastern and central Equatorial Pacific (Trenberth, 1997). In normal, non-El Niño conditions, high CO 2 efflux occurs in the central and eastern Equatorial Pacific due to the upwelling of CO 2 -enriched cold subsurface water by divergence due the trade winds (Feely et al., 2002. In contrast, during El Niños, the easterly trade winds are weakened, leading to decrease in upwelling and an increase of SST. A considerable reduction in outgassing of CO 2 is observed during the El Niño periods (e.g. Feely et al., 2006).
The estimated sea-air CO 2 fluxes in the Equatorial Pacific show a strong correlation with ENSO events. The global anomaly caused by reduced CO 2 efflux to the atmosphere is reinforced by CO 2 uptake in the (Sub)-tropics and the Southern Ocean during El Niño periods (Fig. 4). The combined increase in the extratropical sink (negative flux anomalies) is of similar magnitude to the Equatorial Pacific source decrease. For instance, the average decrease in efflux in the Equatorial Pacific for the El Niño of 1997 and 1998 period (May 1997-May 1998) is 0.17 PgC and the increases in uptake in the (Sub)-tropics and the Southern Ocean are 0.10 and 0.08 PgC, respectively. Intensification of oceanic sink areas outside the Equatorial Pacific during the El Niño periods has also been observed in the sea-air CO 2 fluxes estimated from atmospheric CO 2 inverse modelling (Patra et al., 2005;Rayner et al., 2008). Overall, CO 2 flux variability derived by the ENSO cycle dominates global interannual variability (Fig. 4).

The North Atlantic Oscillation.
The NAO is often defined in terms of a winter index (December-March) of sea level pressure difference between Iceland and Azores because the NAO is most pronounced in amplitude and areal coverage during winter (Hurrell et al., 2003). The NAO exerts a substantial influence on temperature, precipitation, storms and ecosystems of the North Atlantic and surrounding continents. During the positive phase of the NAO, an enhanced gradient in sea level pressure strengthens surface westerly winds over the subpolar gyre. This leads to deeper mixing and a decrease in SST in the subpolar region, and an increase of SST in the western subtropics (Marshall et al., 2001). In the negative and neutral phases, weakened and changed circulations lead to a subpolar warming and subtropical cooling in SST.
For our modelled fluxes from 30 • N to 70 • N in the North Atlantic, there is a strong positive correlation (r = 0.54, P < 0.01) between CO 2 fluxes and NAO index with a 1-year time lag. An even better correlation of 0.70 is shown in Fig. 6 when the fluxes are detrended by removing the mean increase of CO 2 uptake over the three decades. In the positive phase of the NAO, oceanic CO 2 uptake decreases, and in the negative phase of the 1982 1986 1990 1994 1998 2002 2006 Winter NAO index NAO, CO 2 fluxes into the ocean increase in the subtropics and subpolar gyre (Fig. 6).

Year
Our diagnostic model shows predominantly negative pCO 2SW -SST relationships in the subpolar North Atlantic and positive relationships in the subtropics. The SST changes associated with the NAO reinforce the CO 2 fluxes in these gyres and lead to the strong correlation with the NAO index. This positive correlation of CO 2 uptake increase is also shown in an atmospheric CO 2 inversion (Patra et al., 2005) and a biogeochemical general circulation model (Ullman et al., 2009). Ullman et al. (2009) shows an increase of North Atlantic CO 2 sink during the transition period from positive NAO to neutral NAO over mid-1990s to mid-2000s. They attribute this to substantial declines in subpolar gyre convection and vertical inorganic carbon supply related to the phase of the NAO.

The Southern Annular Mode.
The SAM is the principal mode of atmospheric forcing in the Southern Ocean and defined as the difference in mean sea level pressure between 40 • S and 65 • S (Marshall, 2003). The SAM describes the relative atmospheric anomalies at southern mid-and high-latitudes. In positive phase of SAM, sub-polar westerly winds strengthen, leading to surface ocean cooling and increase of upwelling at higher latitudes. Several studies report the connection between the trend in SAM and change in the Southern Ocean sea-air CO 2 flux (Le Quéré et al., 2007;Lovenduski et al., 2008;Metzl, 2009). The Southern Ocean CO 2 sink has been weakened with the recent positive trend of SAM due to increase of outgassing with higher winds in their studies. However, our diagnostic approach shows no significant correlation (P > 0.05) between sea-air CO 2 fluxes and SAM index in the Southern Ocean. This may be due to that variability of pCO 2 in this region is not closely related to SST variation. The validation study of our approach using the output of a GCM also suggests that our empirical approach has limitations in the Southern Ocean.

Assessment of our empirical approach
At local scale, we use time-series measurements made at the Bermuda Atlantic Time Series Study (BATS) and Hawaii Ocean Time Series (HOT) to assess if the subannual pCO 2SW -SST relationships vary from year-to-year at the sites, and to compare fluxes derived from our diagnostic model to those observed at the sites. At regional scale, we compare against a 5-year pCO 2SW record in the North Atlantic from a single SOOP that suggests large interannual changes in the North Atlantic CO 2 flux along the transect . For a global assessment of our method, we apply the empirical method to the SST and pCO 2SW output of a global biogeochemical ocean general circulation model and compare the results with the CO 2 fluxes from the model. The 4 • × 5 • grid cell where BATS is located has a single pCO 2SW -SST relationship for an entire year in our diagnostic model. The derived pCO 2SW -SST relationships from observations at BATS are fairly consistent with our value (Fig. 7a). The 15-year mean of our modelled flux of −1.17 ± 0.16 mol C m −2 yr −1 (1σ ) is also in good agreement with measured sea-air CO 2 fluxes at BATS of −1.18 ± 0.22 mol C m −2 yr −1 (1σ ) (Fig. 8a). Changes in modelled CO 2 fluxes are in phase with those in observations except for the periods of 1995-1997 and 2004-2006, but the magnitude of variability predicted by our empirical method is 30% lower than that of the observations. Two optimum subannual pCO 2SW -SST relationships are derived from pCO 2SW climatology in the grid cell containing the HOT site. One relationship covers wintertime from October to February and the other is from March to September. The yearly March-September relationships determined from measurements at the HOT site are in broad agreement with our relationship (Fig. 7b). However, for the winter season the observed relationships from the HOT data show a large degree of scatter and typically have a larger slope than that estimated from the pCO 2SW climatology. The 16-year average of the CO 2 fluxes at HOT is −0.62 ± 0.20 mol C m −2 yr −1 (1σ ), which is in good agreement with our modelled CO 2 flux of −0.57 ± 0.11 mol C m −2 yr −1 (1σ ) (Fig. 8b). Variations in modelled CO 2 fluxes are in phase with those determined from observations except for the periods of 1998-1999 and 2003-2004. Our modelled variability is significantly lower (45%) than the observed value. Overall, our diagnostic model agrees in magnitude of the average flux with observations but underestimates interannual variability of sea-air CO 2 fluxes in the grid cells including BATS and HOT sites, similar to the results for the previous diagnostic model studies (Lee et al., 1998;Park et al., 2006). Large differences in particular years cause a lower correlation in phase between modelled CO 2 fluxes and observations for the entire period (r = 0.27, P > 0.05). The differences with observations must be due to changes in physical and biological processes that are not strongly correlated with SST anomalies. Substantial interannual variations in nitrogen fixation levels and significant SSS changes from variations of evaporation and precipitation and lateral transport of high salinity water are observed at HOT (Dore et al., 2003(Dore et al., , 2008Keeling et al., 2004). These biological and chemical factors cause variations in pCO 2SW unrelated to SST anomalies that are not accounted for in our model. The difference in scale between the single time-series point and the 4 • × 5 • grid cell probably also plays a role. CO 2 flux variability induced by mesoscale eddies observed at two stations (Bates, 2001; Dore et al., 2008) are not captured by our method. This likely contributes to the lower variability shown by our diagnostic model for the BATS area in spite of good correspondence in the pCO 2SW -SST relationships between model and observations. Thus, the diagnostic approach using 4 • × 5 • grid scale climatological mean values representing a time-space average over the box area does not completely capture the interannual variability in CO 2 fluxes at time-series stations affected by the changes due to local (sub-grid scale) biological and physical processes not closely related to SST. Moreover, sub-grid scale SST variability is not captured by the diagnostic approach. Our approach also does not reproduce the large interannual variations in CO 2 fluxes calculated from SOOP data diagonally crossing the North Atlantic from the Caribbean to the England during the period of . The observations show a change in uptake from 0.39 mol C m −2 yr −1 in 2002 to a maximum uptake of 0.95 mol C m −2 yr −1 in 2005 and a decreasing trend from then on. Our results for the overlapping grid cells show the same trend from 2003 to 2007 as the observations but with a much smaller magnitude (Fig. 9).
3.5.2. Validation of our approach using a global biogeochemical ocean model. Because there are only limited timeseries observations to validate our approach, the veracity of this empirical model was tested using the output of a 3-D biogeochemical ocean general circulation model (GCM). Model outputs including pCO 2 , pCO 2SW and SST are obtained from a multi-decade hindcast simulation   biological-chemical module, the Biogeochemical Elemental Cycle model (Doney et al., 2009a,b). The model output is regridded to the same spatial grid (4 • × 5 • ) as the global pCO 2 climatology. In this GCM, pCO 2SW values are computed from the model prognostic variables, total dissolved inorganic carbon, total alkalinity, SST and SSS, using a full seawater carbonate thermodynamic code. Subannual pCO 2SW -SST relationships are derived from year 2000 model pCO 2SW and SST data in the same way as was done using the pCO 2 climatology of Takahashi et al. (2009a). The mean correlation coefficient and standard deviation of the pCO 2SW -SST relationships derived from the GCM output for all the grid cells is 0.89 ± 0.13, similar to that of 0.83 ± 0.14 of our empirical approach using the pCO 2 climatology. Global net sea-air CO 2 fluxes for the period of 1982-2006 are estimated from these subannual relationships with model SST anomalies and monthly NCEP-II wind speed data. In the central and eastern Equatorial Pacific (6 • N-10 • S and 80 • W-165 • E), the same empirical pCO 2SW -SST equations derived from in situ measurements are used as in our empirical approach using observations. The global sea-air CO 2 fluxes are also calculated directly from the pCO 2 produced by the GCM and monthly NCEP-II wind speeds. Figure 10a shows the comparison of anomalies in global net sea-air CO 2 fluxes between the GCM output (PROG GCM) and our diagnostic approach applied to the GCM output (DIAG GCM). They show overall good agreement in terms of phase, but there are significant differences in global CO 2 fluxes for several years (1986-1987, 1997, 2003-2006). Except for those years, our diagnostic approach applied to GCM output reproduces interannual variability in global sea-air CO 2 fluxes with similar magnitude and sign as the prognostic GCM results (Fig. 10b). The 25-year global mean sea-air CO 2 flux and its interannual 1982 1986 1990 1994 1998 2002 2006 CO 2 flux anomaly (Pg C yr  (Takahashi et al., 2009a). This is the same as the black line in Fig. 4, except that the time series is shorter here. (b) Correlation of global CO 2 flux anomalies between the DIAG GCM and PROG GCM. The dotted line shows the regression with slope of 1. Open circles indicate the periods from 1986-1987, 1997 and 2003-2006 when CO 2 flux anomalies show large differences between two results in (a) and these circles are excluded in the regression (solid line).

Year
variability are −1.33 ± 0.17 PgC yr −1 (1σ ) for the prognostic GCM output compared to −1.28 ± 0.14 PgC yr −1 (1σ ) for our diagnostic approach applied to the GCM output ( Table 2). The 23% lower interannual variability using the empirical approach applied to the GCM output is attributed to a few specific regions (Fig. 11). The largest difference is found in the Southern Ocean. The prognostic GCM output shows large interannual variability in the latitude band between 40 • S and 60 • S (Fig. 11). Our estimates using optimum subannual pCO 2SW -SST relationships obtained from GCM output and SST anomalies accounts for only half of interannual variability of model output in the Southern Ocean (±0.09 PgC yr −1 ) ( Fig. 12d; Table 2). Many grid cells in the Southern Ocean show low correlations (P > 0.05) in annual The percentage of CO 2 variability relative to the regional 25-year mean annual CO Fig. 11. Global distributions of interannual variability magnitude estimated from (a) the GCM output (PROG GCM) and (b) our diagnostic approach applied to the GCM pCO 2SW data (DIAG GCM). The magnitude of interannual variability is expressed as a standard deviation (1σ ) of the annual net sea-air CO 2 fluxes in each grid cell.  Table 2. net sea-air CO 2 fluxes between our diagnostic approach applied to the GCM output and the prognostic GCM output (Fig. S4) indicating that pCO 2 variability in the Southern Ocean in the model is driven by factors not directly related to SST.
Sea-air CO 2 fluxes estimated from the empirical pCO 2SW -SST equations applied to GCM SST in the central and eastern Equatorial Pacific agree well with prognostic GCM outputs in terms of both phase and magnitude (Figs 12a and S4) as was noted in a similar comparison by Doney et al. (2009a). The results obtained from our approach in northern high-latitude oceans also show overall good agreement with GCM output fluxes ( Fig. 12c; Table 2). However, there are differences in the patterns of variability at the sub-basin scale. For the prognostic model output, the western North Pacific shows larger interannual variability than in the surrounding areas, while the eastern North Pacific has larger variability for our diagnostic approach applied to the GCM output (Fig. 11). At high latitudes in the North Atlantic, our diagnostic approach applied to the GCM output shows large interannual variability over the entire region but the prognostic GCM output only shows large variability in the western part (Fig. 11). The variability in annual net sea-air CO 2 fluxes reproduced by our diagnostic approach are generally in phase with the prognostic GCM output in the Atlantic Ocean (Fig. S4). In the subtropical and tropical oceans, our approach applied to the GCM output has about 30% higher magnitude of variability in CO 2 fluxes compared to the prognostic GCM output. There are large differences between the two approaches for the period of 2003-2006 (Fig. 10a). The difference during this time is mainly caused by smaller CO 2 uptake in the subtropical and Southern oceans by our diagnostic approach (Figs 12b and  d). During the El Niño periods of 1986-1987 and 1997, our diagnostic approach applied to the GCM output exhibits much larger negative global flux anomalies (net ocean uptake) than that for the prognostic GCM output (Fig. 10a). The different model CO 2 flux anomaly curves are quite similar in the Equatorial Pacific (Fig. 12a), and the discrepancy between the prognostic GCM output and the diagnostic approach during these time periods arises in the extratropics, in particular the Southern Ocean (Fig.  12d). These differences for the specific periods contribute to the lower correlation in global net sea-air CO 2 fluxes as well as in regional fluxes of the subtropical and Southern oceans ( Table 2).
The GCM flux output (PROG GCM) and our diagnostic approach applied to the GCM output (DIAG GCM) are compared with the global net sea-air CO 2 flux anomalies based on pCO 2 climatology and our diagnostic approach (DIAG CLIM) in Fig.  10a. The approaches generally agree in terms of phase, except for the period of 1987-1989. During this transition period from El Niño to La Niña, when global CO 2 fluxes into the ocean decrease in DIAG GCM and DIAG CLIM while the PROG GCM shows nearly constant positive CO 2 flux anomaly (Fig. 10a).
Significant differences are observed when comparing the influence of NAO variability on the CO 2 fluxes from the GCM with that of our approach using pCO 2 climatology. In an analysis of the GCM output, Thomas et al. (2008) found a significant correlation of sea-air CO 2 flux with the NAO index only for the western subpolar gyre. However, their correlation shows the opposite phase with our result using pCO 2 climatology and SST anomalies. In the positive phase of the NAO index, the GCM shows increase in CO 2 uptake while our approach shows the opposite (Fig. 6). This is because our diagnostic model only accounts for changes related to SST while the GCM includes other processes controlling pCO 2SW that are not related to SST variations (Doney et al., 2009a). Thomas et al. (2008) explained this correlation with North Atlantic Current (NAC) transporting water with low CO 2 concentration into the region. Although several studies suggest that variability of oceanic CO 2 uptake in the North Atlantic is related to the NAO, the phase, location and magnitude and even direction of CO 2 fluxes related with this climate mode differ greatly (Corbière et al., 2007;Schuster and Watson, 2007;Thomas et al., 2008;Schuster et al., 2009;Ullman et al. 2009). Continued basin-wide observations are warranted to resolve these discrepancies.
The mechanisms driving the interannual variability in net sea-air CO 2 flux simulated by the GCM vary by region (Doney et al., 2009a). Changes in circulation are a dominant control on variability in the Equatorial Pacific, high northern oceans and Southern Ocean. In subtropics and tropical Atlantic, physical chemical processes are dominant. The variation of dust deposition is also key factor causing large interannual variability in CO 2 fluxes in the western North Pacific and the Southern Ocean in the model, accounting for up to 50% of total regional CO 2 flux variability. Iron supply to high nitrate and low chlorophyll (HNLC) regions through atmospheric dust deposition downwind of desert source regions encourages biological production, which in turn reduces pCO 2SW and increases net oceanic CO 2 uptake. Enhanced productivity by dust deposition will not be captured by our empirical approach. The ability of our approach to capture much of the variability shown in the GCM suggests that it has appreciable skill on global scale. However, there are also differences that show up when applying our method to a GCM, clearly showing that all variations in biogeochemistry controlling pCO 2 cannot be captured just using SST anomalies.

Conclusions
The interannual variability in net sea-air CO 2 fluxes is estimated from subannual pCO 2SW -SST relationships approximated by segmented linear fits to new monthly climatological maps of pCO 2SW . The optimum subannual regressions are robust over most of the ocean with a global average correlation coefficient (R 2 ) for all grid cells of 0.83. In the comparison with observations, our approach does not fully capture the magnitude of interannual variability in CO 2 fluxes due to local scale drivers unrelated to SST changes, and thus, may underestimate the variability. A comprehensive evaluation is performed by applying it to output of global biogeochemical ocean general circulation model. The comparison shows that our approach using subannual pCO 2SW -SST relationships and interannual SST anomalies to predict interannual changes in pCO 2SW reasonably captures the interannual variability of CO 2 fluxes in the GCM in most ocean regions except the Southern Ocean and for select time periods. Our diagnostic approach captures the phase and magnitude of variability for most regions and time periods albeit with an about 23% lower amplitude when the derived pCO 2SW -SST relationships fairly represent actual mean state. The CCMP wind product has 20% greater variability than the smoothed NCEP-II wind product, which shows that the wind speed product can be equally important as accurate estimate of pCO 2 in the study of CO 2 flux variability.
The estimated oceanic CO 2 variability is closely related to prominent climate modes such as the NAO and the ENSO. The good predictability of the variation in fluxes due to ENSO is in large part due to sustained observations of pCO 2SW (and SST) in the region and the adaptive relationships between pCO 2SW and SST that change with phase of ENSO and over time . Our empirical approach can be improved with increasing high quality pCO 2SW observations to create regional regressions of pCO 2SW and SST, much like is done for the Equatorial Pacific rather than relying on the climatology that invariably has smoothed SST and pCO 2SW seasonal variations. Moreover, sustained long-term observations covering large areas will allow us to better assess the validity of our diagnostic at regional scales.

Acknowledgments
This work was supported by the Korea Research Foundation Grant funded by the Korean Government (KRF-2008-357-C00167), the NOAA Global Carbon Cycle Program (GC07-193 and NA07OAR4310098) and the NOAA Climate Observations Division. This research was carried out, in part, under the auspices of the Cooperative Institute for Marine and Atmospheric Studies (CIMAS), a Joint Institute of the University of Miami and the National Oceanic and Atmospheric Administration (cooperative agreement No. NA17RJ1226). We are grateful to Nick Bates, John Dore and other scientists at the BATS and HOT sites for contributing their time-series data.

Supporting information
Additional supporting information may be found in the online version of this article: Table S1. Comparison of mean net sea-air CO 2 fluxes and variability (1σ ) estimated from three fixed seasonal and optimum subannual relationships. Fig. S1. Map of monthly distribution of pCO 2SW -SST relationships for each grid cell. Fig. S2. Map of monthly distribution of R 2 of pCO 2SW -SST relationships for each grid cell. Fig. S3. Difference in interannual variability estimated from three fixed seasonal and optimum subannual relationships. Fig. S4. Correlation of annual net sea-air CO 2 fluxes between DIAG GCM and PROG GCM.
Please note: Wiley-Blackwell is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.