# Global surface-ocean pCO2 and sea-air CO2 flux variability from an observation-driven ocean mixed-layer scheme

Christian Rödenbeck, Ralph F. Keeling, Dorothee C. E. Bakker, Nicolas

Metzl, Are Olsen, Christopher L. Sabine, M. Heimann

► To cite this version:

Christian Rödenbeck, Ralph F. Keeling, Dorothee C. E. Bakker, Nicolas Metzl, Are Olsen, et al. Global surface-ocean pCO2 and sea-air CO2 flux variability from an observation-driven ocean mixed-layer scheme. Ocean Science, 2013, 9, pp.193-216. 10.5194/OS-9-193-2013. hal-00810525

HAL Id: hal-00810525

https://hal.science/hal-00810525

Submitted on 4 Dec 2015

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# Global surface-ocean $p^{CO_2}$ and sea-air $CO_2$ flux variability from an observation-driven ocean mixed-layer scheme

C. Rödenbeck¹, R. F. Keeling², D. C. E. Bakker³, N. Metzl⁴, A. Olsen⁵,⁶,⁷, C. Sabine⁸, and M. Heimann¹

¹Max Planck Institute for Biogeochemistry, Jena, Germany

²Scripps Institution of Oceanography, University of California, San Diego, USA

³School of Environmental Sciences, University of East Anglia, Norwich Research Park, Norwich, UK

⁴LOCEAN-IPSL, CNRS, Paris, France

⁵Institute of Marine Research, Bergen, Norway

⁶Uni Bjerknes Centre, Bergen, Norway

⁷Bjerknes Centre for Climate Research, Bergen, Norway

⁸NOAA Pacific Marine Environmental Laboratory, Seattle, USA

*Correspondence to:* C. Rödenbeck (christian.roedenbeck@bgc-jena.mpg.de)

Received: 27 April 2012 – Published in Ocean Sci. Discuss.: 27 June 2012

Revised: 23 January 2013 – Accepted: 29 January 2013 – Published: 1 March 2013

**Abstract.** A temporally and spatially resolved estimate of the global surface-ocean $CO_2$ partial pressure field and the sea-air $CO_2$ flux is presented, obtained by fitting a simple data-driven diagnostic model of ocean mixed-layer biogeochemistry to surface-ocean $CO_2$ partial pressure data from the SOCAT v1.5 database. Results include seasonal, interannual, and short-term (daily) variations. In most regions, estimated seasonality is well constrained from the data, and compares well to the widely used monthly climatology by Takahashi et al. (2009). Comparison to independent data tentatively supports the slightly higher seasonal variations in our estimates in some areas. We also fitted the diagnostic model to atmospheric $CO_2$ data. The results of this are less robust, but in those areas where atmospheric signals are not strongly influenced by land flux variability, their seasonality is nevertheless consistent with the results based on surface-ocean data. From a comparison with an independent seasonal climatology of surface-ocean nutrient concentration, the diagnostic model is shown to capture relevant surface-ocean biogeochemical processes reasonably well. Estimated interannual variations will be presented and discussed in a companion paper.

## 1 Introduction

The oceans are considered the dominant player in the global carbon cycle on long timescales, e.g. in the glacial-interglacial cycles (e.g. Sigman and Boyle, 2000). On a multi-millennial timescale, the oceans will be the sink for 80–95 % of the anthropogenic $CO_2$ emissions, and 70–80 % on a timescale of several hundred years (Archer et al., 1997). Currently, the oceans take up about 25 % of the emissions (Sarmiento et al., 2010). Concerns exist, however, that the sink efficiency may decrease in the coming decades as a consequence of anthropogenic climate change, as suggested by model projections (e.g. Sarmiento and Le Quéré, 1996; Matear and Hirst, 1999; Joos et al., 1999) and tentatively confirmed by data analysis (e.g. Le Quéré et al., 2007, 2010). As a prerequisite to understanding the involved processes, one needs to quantify sea-air $CO_2$ fluxes, their variability, and their response to forcing.

Currently, two data streams are used to estimate the variability of sea-air $CO_2$ fluxes:

- Based on measurements of the $CO_2$ partial pressure ($p^{CO_2}$) of surface water, the sea-air $CO_2$ flux is calculated through a gas exchange parameterization. As $p^{CO_2}$ measurements only exist at discrete locations and ship tracks, spatial and temporal interpolation is needed. The following interpolation methods have been
proposed: (1) statistical interpolation combined with an advection–diffusion equation (Takahashi et al., 2009); (2) purely statistical interpolation with error quantification (Jones et al., 2012b); (3) multi-linear regressions between $p^{CO_2}$ and ocean state variables (e.g. Park et al., 2010; Chen et al., 2011); (4) neural networks learning the relationship between $p^{CO_2}$ and oceanic state variables available from remote sensing platforms and ocean reanalysis projects (e.g. Lefèvre et al., 2005; Telszewski et al., 2009); (5) assimilation of the $p^{CO_2}$ data into a process model of ocean biogeochemistry (Valsala and Maksyutov, 2010; While et al., 2012). These methods are complementary in certain aspects: the more statistical methods are dominated by the data themselves but are strongly affected by spatial and temporal gaps, while the more complex methods more easily spread the information but are dependent on driver data sets or even on model formulations.

- The CO₂ exchange between the Earth surface and the atmosphere can be quantified based on atmospheric CO₂ mixing ratio measurements (e.g. Conway et al., 1994) by the atmospheric transport inversion (e.g. Newsam and Enting, 1988; Rayner et al., 1999; Bousquet et al., 2000; Rödenbeck et al., 2003; Baker et al., 2006) (see Sect. 2.1). This estimation does not involve any parameterization of gas exchange. However, ocean fluxes are difficult to detect by this method in most parts of the globe, because their imprint on the atmospheric mixing ratio records is small compared to that of the much more variable land fluxes: Even if ocean-internal processes (biology, transport) cause mixed-layer carbon sources and sinks of comparable variability as the land biosphere, the resulting sea-air CO₂ exchange is much less variable and smoothed out in time because the carbonate chemistry of seawater slows the equilibration rate of dissolved CO₂ with the atmosphere. In addition, the responses to warming/cooling partially counteract the effects of ocean-internal sources/sinks. A further problem of atmospheric inversions based on data from a discrete set of measurement sites is that information is only provided on scales comparable to or larger than the distance between the sites or the sampling frequency, while variability exists also on smaller scales.

A further data-based method to estimate sea-air CO₂ fluxes uses ocean-interior carbon data in inverse calculations of oceanic transport (Gloor et al., 2003; Mikaloff Fletcher et al., 2006). This method is independent of parameterizations of gas exchange as well. However, it only yields long-term sea-air CO₂ fluxes over large spatial regions, not its temporal or high-resolution spatial variability. The long-term global sea-air CO₂ flux has also been estimated from observed trends in atmospheric oxygen (Keeling and Shertz, 1992) as well as from ¹³C isotopic ratios in atmospheric CO₂ (Ciais et al., 1995).

In view of the above-mentioned problem of the atmospheric CO₂ inversions, most of these studies use the sea-air flux climatology by Takahashi et al. (2009) (or earlier versions) as a Bayesian prior (in some cases, the long-term fluxes from ocean-interior inversions are used in the prior as well, e.g. Rödenbeck et al. (2003), or formalized as joint inversion by Jacobson et al. (2007)). Out of this context, the database sea-air CO₂ fluxes presented here have been developed aiming to (1) provide not only a monthly seasonal cycle but variability also on short-term (daily) and interannual timescales, (2) use an assimilation scheme that can easily and self-consistently compare or even combine the observational information of atmospheric and surface-ocean data, and (3) use a framework that can be extended to also incorporate constraints from other tracers, such as oxygen.

The presented results are based on the newly available Surface Ocean CO₂ Atlas (SOCAT) database (version 1.5) of CO₂ fugacity measurements (Pfeil et al., 2012, http://www.socat.info/). We describe the estimation method and test its performance, also using independent data. As a first step of analysis, this study mainly focuses on the mean seasonal cycle, because on this timescale (1) we can compare to the widely used monthly $p^{CO_2}$ climatology by Takahashi et al. (2009), (2) we can investigate mutual consistency between the surface-ocean $p^{CO_2}$ constraint and that by atmospheric CO₂ data as the signal/noise ratio is large, and (3) we can test the plausibility of our process representations invoking a seasonal surface-ocean phosphate (PO₄) climatology. Estimated interannual variations of the sea-air CO₂ flux are presented and discussed in a companion paper (Rödenbeck et al., 2013).

## 2 Method

### 2.1 Concept – overview

In a classical atmospheric inversion, a spatio-temporal field of surface-to-atmosphere carbon fluxes is estimated such that its corresponding mixing ratio field – as simulated by an atmospheric tracer transport model – matches as closely as possible a set of mixing ratio observations. The match is gauged by a quadratic cost function to be minimized (Appendix A2).

Here we extend the inversion framework by not only considering the process of atmospheric transport but also processes in the oceanic mixed layer: The atmospheric transport model is supplemented by a chain of parameterizations of gas exchange, carbonate chemistry, and a carbon budget equation, which determine the sea-to-air CO₂ exchange ($f_{ma}^{CO_2}$) as a function of ocean-internal carbon sources and sinks (Fig. 1). Then the inversion is not directly adjusting the sea-to-air fluxes any more, but only indirectly through adjusting the ocean-internal fluxes instead (Fig. 2, from left to middle). This offers two advantages: (1) the equations describing
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Summary of the main quantities and process parameterizations involved in the diagnostic ocean mixed-layer model. Thick boxes denote the process parameterizations given in Sect. 2.2, each expressing the quantity right above as a function of the quantity right below. Quantities at the arrows on the left represent driver data entering the parameterizations. See Table 2 for mathematical symbols, and Appendix A1 for the equations and further explanation.</figcaption></figure>

the individual processes impose spatial and temporal structure to the sea-air fluxes, based on the spatio-temporal information from the driver data (sea-surface temperature (SST), wind speed, etc., Fig. 1). (2) The explicit representation of oceanic processes involves further quantities, notably sea-surface CO₂ partial pressure. Through cost function contributions gauging the match of modelled and measured pCO₂, these data can be used as an observational constraint replacing the atmospheric data (Fig. 2, right). In this mode of operation, atmospheric transport models and atmospheric data are actually no longer used in the inverse calculation, but the Bayesian framework, including a-priori spatial and temporal correlations, is still applied in the same way as in the atmospheric mode.

The calculation is global over the time period 1985–2011, with a temporal resolution of 1 day and a horizontal resolution of approximately 4° latitude × 5° longitude (TM3 model grid).

## 2.2 Process parameterizations

All processes considered explicitly are given in Fig. 1 and summarized in the following. Details, including the equations used, are found in Appendix A.

**Atmospheric transport.** Atmospheric CO₂ mixing ratio fields in response to surface-to-atmosphere fluxes are simulated by the global off-line atmospheric transport model TM3 (Heimann and Körner, 2003) with a spatial resolution of ≈ 4° lat. × 5° long. × 19 vertical levels. The model is driven by 6-h interannual meteorological fields derived from the National Centers for Environmental Prediction (NCEP) reanal-

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Illustration of the inverse procedure (three modi). Boxes denote the process parameterizations from Fig. 1, causally linking quantities from bottom to top. Double arrows symbolize the matching between observed and modelled quantities, as gauged by a cost function *J* (Appendix A2). Thin arrows indicate the adjustments of the respective unknowns done to minimize the model-data mismatch. Left: pure atmospheric transport inversion not used here but given for reference: the sea-air fluxes *f*<sub>ma</sub><sup>CO2</sup> are adjusted to match atmospheric mixing ratios (terrestrial net ecosystem exchange is adjusted as well but omitted from this schematic for clarity). Middle: in the case of **ATM**, ocean-internal sources/sinks *f*<sub>int</sub><sup>DIC</sup> are adjusted instead of sea-air fluxes to match atmospheric CO₂ mixing ratios (again net ecosystem exchange is further adjusted). Right: in the case of **SFC**, ocean-internal sources/sinks are adjusted to match the surface-ocean CO₂ partial pressure observations. Sea-air fluxes are then calculated from the estimated partial pressure field.</figcaption></figure>

ysis (Kalnay et al., 1996). The model fields are sampled at the location and time of the individual mixing ratio measurements used. TM3 has performed well in intercomparisons of state-of-the-art global tracer transport models (e.g. Stephens et al., 2007; Law et al., 2008).

**Solubility and gas exchange.** Diffusive sea-to-air gas exchange is proportional to the over-/undersaturation of CO₂ in the surface ocean and to the piston velocity with the quadratic wind speed dependence by Wanninkhof (1992), scaled to match the global average piston velocity of Naegler (2009) (Appendix A1.1).

**Carbonate chemistry.** The carbon species relevant for gas exchange (CO₂) only account for a small part of the carbon relevant in the ocean-internal budget (dissolved inorganic carbon, DIC). The link between CO₂ abundance (expressed in terms of partial pressure *p*<sub>m</sub><sup>CO2</sup>) and DIC abundance (in terms of its concentration *C*<sub>m</sub><sup>DIC</sup>) is determined by chemical equilibria, which we assume to be attained instantaneously. The non-linear dependence of *p*<sub>m</sub><sup>CO2</sup> on *C*<sub>m</sub><sup>DIC</sup> has been linearized in standard ways (Sarmiento and Gruber, 2006). The chemistry parameterization also contains a temperature-dependent factor, and contributions from seasonal variations in alkalinity and salinity (Appendix A1.2).
*Mixed-layer carbon budget.* Changes in the spatio-temporal field of dissolved inorganic carbon in the ocean mixed layer need to be balanced by the sum of fluxes (Appendix A1.3). As the sea-air flux itself depends on carbon concentration, this balance can be expressed as a linear first-order differential equation. Its solution gives DIC concentration $C_m^{DIC}$ as a function of the ocean-internal sinks $f_{int}^{DIC}$ (representing the total effect of biological conversion and vertical/horizontal advection/diffusion on mixed-layer carbon concentration). An additional important item in the budget is the re-entrainment during mixed-layer deepening of carbon left behind previously during mixed-layer shoaling, which we represent as a “history flux” $f_{hist}$. We further consider the influence of freshwater fluxes.

## 2.3 Main adjustable degrees of freedoms

The ocean-internal carbon sources and sinks ($f_{int}^{DIC}$) at the end of the described chain of process parameterizations is the basic unknown to be adjusted to match the surface-ocean $p^{CO_2}$ or atmospheric $CO_2$ data (next section). It is treated similarly to the unknown sea-air flux in the pure atmospheric transport inversion (Fig. 2), including Bayesian a-priori spatial and temporal correlations. The detailed specification and further explanations are found in Appendix A2.2.

## 2.4 Data constraints – base runs

We will present results of two cases that differ in the data set used as main constraint (Fig. 2):

**SFC.** In the main case used to create the primary product of this study, the ocean part of the diagnostic scheme is fit to surface $p^{CO_2}$ data points from the SOCAT v1.5 database (Pfeil et al., 2012, http://www.socat.info/). Data pre-treatment and further details are given in Table 1. SOCAT data cover about 6 million pixels/time steps globally within the calculation period (until 2007, see Supplement Figs. S7.3 and S7.4 for data distribution).

**ATM.** As a comparison case used to assess the consistency between surface-ocean and atmospheric data, the diagnostic scheme is fit to atmospheric $CO_2$ mixing ratios measured approximately weekly or hourly at a set of sites by various institutions. Details are given in Table 1.

# 3 Results and discussion

The main product of this study is a data-based estimate of the global spatio-temporal $CO_2$ partial pressure field and ultimately the sea-air $CO_2$ fluxes. Here we characterize and evaluate its robustness, errors, and information content. In addition, we compare these results with independent data and previous estimates. The consideration of interannual variations is done in the companion paper (Rödenbeck et al., 2013).

## 3.1 Overview

To illustrate the characteristics of the quantities and the temporal scales involved in the scheme, Fig. 3 shows time series of key quantities (as of run SFC, blue) at the example pixel around Station M in the Norwegian Sea (66°N, 2°E). The ocean-internal sources/sinks (bottom panel) are relatively smooth, as the a-priori temporal correlations (Appendix A2.2) suppress fast variations. The mixed-layer DIC concentration (panel above) responds to these ocean-internal fluxes and the emerging sea-air exchange. Being their temporal integral, the DIC concentration is slightly shifted in phase with respect to the ocean-internal fluxes. In addition, the DIC concentration is rising in response to the atmospheric $CO_2$ increase. The surface-ocean $CO_2$ partial pressure (3rd panel from bottom) shares the rise and the variations with the DIC concentration, but the seasonality is again slightly shifted because of the temperature and alkalinity effects. Finally, the sea-air $CO_2$ flux (top panel) is dominated by the $p^{CO_2}$ variability; in addition, it shows high-frequency (daily) variability due to wind speed and solubility changes. Figure 3 also illustrates the dampening effect of the carbonate chemistry on sea-air exchange, as the seasonal amplitude of the sea-air flux (top) is much lower than that of the ocean-interior flux (bottom).

The spatial resolution and domain of the calculation is illustrated by Fig. 4, showing the amplitude of the mean seasonal cycle of $p^{CO_2}$ for each pixel.

## 3.2 Data and model constraints

The estimated $p^{CO_2}$ and sea-air $CO_2$ flux fields combine information from the partial pressure data and from the process parameterizations. To illustrate this, Fig. 3 also shows the a-priori state (thin dashed grey) defined by $f_{int, pri}^{DIC} = 0$ (Appendix A2.2, see bottom panel). The choice $f_{int, pri}^{DIC} = 0$ means that, without data knowledge, it is equally likely to have an internal source or an internal sink at any given location and time, and thus corresponds to a state of no information on $f_{int}^{DIC}$. The variations in the prior values of $C_m^{DIC}$, $p^{CO_2}$, and the sea-air fluxes $f_{ma}^{CO_2}$, which then follow from the process parameterizations, only comprise responses to variations in sea surface temperature and the other driving variables (including rising atmospheric $CO_2$), but miss variations in the ocean-internal sources/sinks. As this a-priori $p^{CO_2}$ field contradicts the data (black dots), the estimation procedure now adjusts the internal flux $f_{int}^{DIC}$ in such a way that the data are matched as closely as possible (blue). Most prominently, this reduces or – as at the example pixel of Fig. 3 – even reverses the seasonal variations of $p^{CO_2}$ (reflecting that thermal, biological, and physical effects partially oppose each other). The fit at further locations is given below (Sect. 3.4).

For most pixels, the SOCAT data set only sporadically contains data points (in particular, periods of dense
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Example time series (full daily resolution) at the model pixel around Station M (66° N, 2° E) of the main quantities of the diagnostic model (Fig. 1): sea-air $CO_2$ fluxes ($f_{ma}^{CO_2}$), surface-ocean $CO_2$ partial pressure ($p_m^{CO_2}$), mixed-layer DIC concentration ($C_m^{DIC}$), as well as ocean-internal carbon addition to (removal from) the mixed layer ($f_{int}^{DIC}$). Estimates have been obtained by fitting the diagnostic scheme to SOCAT surface-ocean $CO_2$ partial pressure data (SFC, blue). The SOCAT data points that happen to fall in the example pixel are shown as black dots in the $p_m^{CO_2}$ panel. The Bayesian prior ($f_{int}^{DIC} = 0$ and corresponding $C_m^{DIC}$, $p_m^{CO_2}$, and $f_{ma}^{CO_2}$, Appendix A2.2, Sect. 3.2) is also given (thin dashed grey lines).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Amplitude of the seasonal cycle of surface-ocean $CO_2$ partial pressure (µatm) estimated by fitting the diagnostic scheme to the SOCAT data (run SFC). The amplitude is given as temporal standard deviation of the 1997–2009 monthly mean $p_m^{CO_2}$ at each pixel.</figcaption></figure>
<figure><img src="image_8.png" /></figure>

**Fig. 5.** Illustration of the information flow in the spatio-temporal extrapolation through a-priori correlations. Where data values exist, pixels are directly constrained (left column, which is identical to Fig. 2). Due to the a-priori spatial correlations (Appendix A2.2), adjustments to $f_{\text{int}}^{\text{DIC}}$ at the constrained pixel also change $f_{\text{int}}^{\text{DIC}}$ at neighbouring unconstrained pixels. $p^{\text{CO}_2}$ is then calculated passively from the extrapolated $f_{\text{int}}^{\text{DIC}}$ by the process parameterizations, using the local values of the driving fields (right column). In data-dense areas, local adjustments to $f_{\text{int}}^{\text{DIC}}$ (both for constrained and for extrapolated pixels) will need to compromise between several data constraints within the correlation radius (de-weighted with distance due to the decay of the correlation, Fig. A1). As with any smoothing, this may be beneficial in suppressing outliers, but also adverse in partly suppressing small-scale signals (the structures from the driving data of the parameterizations are never smoothed however, as the smoothing only acts on $f_{\text{int}}^{\text{DIC}}$). In data-poor areas, there may exist pixels that do not have any constrained pixel within the correlation radius; then $f_{\text{int}}^{\text{DIC}}$ will stay close to the prior. In time, extrapolation happens not only due to the a-priori correlations, but additionally due to the memory effect resulting from the relaxation time of the budget equation Eq. (A21).

regular sampling, as at Station M since 2006, are a rare exception). Thus, only a small fraction of pixels/time steps is constrained directly in the described way. However, most parts of the space-time domain are still constrained indirectly through extrapolation via the spatial and temporal a-priori correlations defined in Appendix A2.2. The mechanism is illustrated in Fig. 5. In particular, as seasonal variations in $f_{\text{int}}^{\text{DIC}}$ are allowed to be adjusted much more freely than non-seasonal variations (equivalent to a-priori correlations between adjustments to $f_{\text{int}}^{\text{DIC}}$ at any given day of year in each year, Appendix A2.2), the mean seasonal cycle is extrapolated throughout the time period of the calculation. This also means that data points in any year can contribute to constrain this mean seasonality. Thus, though there are regions or periods too far away from data points for extrapolation, larger-scale seasonal variations are constrained almost everywhere in the ocean. This has been confirmed by assessing how well a given $p^{\text{CO}_2}$ field can be retrieved by the scheme

from pseudo data sampled from this field at the times and locations of the SOCAT data (Appendix B). Constraining inter-annual variations requires more even data distribution in time than constraining seasonal variability; this density of data is only available in some areas of the ocean (Rödenbeck et al., 2013).

The spatio-temporal a-priori correlations also act to smooth $f_{\text{int}}^{\text{DIC}}$ on small scales (Fig. 3, bottom), as they damp small-scale and short-term variations. This has been done because the available data do not have the spatial or temporal density required to constrain these fine-scale features adequately. Nevertheless, despite the smooth $f_{\text{int}}^{\text{DIC}}$ the corresponding $p^{\text{CO}_2}$ and sea-air $CO_2$ flux fields (Fig. 3, top) do comprise fast variations as represented in the process parameterizations, including responses to variations in temperature (changes in solubility and chemical equilibrium) or in wind speed (e.g. accelerated depletion of an oversaturated mixed layer in a high wind speed event, Bates et al., 1997). These effects already account for a considerable part of short-term variations (see Sect. 3.4 below). Only variations related to fast ocean-internal processes, such as the sub-weekly variations of upwelling events or algal blooms, will be missing. In any case, the pseudo-data tests (Appendix B) confirm that the degrees of freedom in the diagnostic scheme provide sufficient flexibility to follow variability as contained in the $p^{\text{CO}_2}$ fields from a biogeochemical process model simulation (as the model has less day-to-day variability than the real $p^{\text{CO}_2}$ field, however, this test may underestimate errors from aliasing such high-frequency variations into variations on the resolved timescales).

### 3.3 Robustness

Figures 6 and 7 show the seasonality of the sea-air $CO_2$ flux and the surface-ocean $CO_2$ partial pressure as estimated from SOCAT data by run SFC. The fields have been integrated/averaged over a set of regions splitting the ocean into basins and latitude bands.

To test how robustly the results are determined from the data and to assess the impact of errors in the parameterizations and their driving fields, individual important model elements have been changed in sensitivity runs: (i) increase and decrease of the a-priori uncertainty by a factor 2, leaving more/less freedom for inverse adjustments, (ii) decrease of the a-priori uncertainty of non-seasonal variations in $f_{\text{int}}^{\text{DIC}}$ by a factor of 2, (iii) increase in the spatial correlation lengths by a factor of 3, (iv) increase and decrease of the global mean piston velocity by 3.2 cm (range given by Naegler, 2009), and (v) increase and decrease mixed-layer depth $h$ by a factor of 2. The resulting range of results (grey bands in Figs. 6 and 7) is generally very narrow compared to the seasonal amplitude, in particular for $p^{\text{CO}_2}$. In most regions, the sensitivity of the sea-air flux is slightly larger than that of $p^{\text{CO}_2}$, almost entirely due to the uncertainty of gas exchange (variation in piston velocity). Larger spread in $p^{\text{CO}_2}$ than in the sea-air flux
<figure><img src="image_9.png" /></figure>

**Fig. 6.** Monthly mean seasonal cycle of sea–air $CO_2$ fluxes, estimated by fitting the diagnostic scheme to SOCAT data (run SFC, blue). The grey shading around the standard result comprises sensitivity cases listed in Sect. 3.3. Fluxes have been integrated over latitude bands split into ocean basins (including open ocean, coastal areas and marginal seas; panels in geographical arrangement). Seasonal cycles have been calculated over the period 1997–2009; the first half year is repeated for clarity.

is found mostly in the North Pacific and the temperate North Atlantic, entirely due to a few marginal areas (including pixels in the far north of the Pacific, or the Black Sea and the Caspian Sea added to the Atlantic region, see Fig. 4) where the $p^{CO_2}$ seasonality is exceptionally high due to high SST seasonality, but where gas exchange is low such that their contribution to regional sea–air flux is nevertheless small. If $p^{CO_2}$ is averaged only over the open ocean, the $p^{CO_2}$ spread becomes narrow in all regions (see Supplement Fig. S7.8).

## 3.4 Fit to the data

The match of the estimated $p^{CO_2}$ field and measured data points somewhat depends on the considered location (Fig. 8). At the North Pacific example pixel (i.e. within the region of largest spread among the sensitivity cases in Fig. 7), the seasonality is underestimated by some µatm. This may result

<figure><img src="image_10.png" /></figure>

**Fig. 7.** Analogously to Fig. 6, monthly mean seasonal cycle of surface-ocean carbon dioxide partial pressure $p_m^{CO_2}$, estimated by fitting the diagnostic scheme to SOCAT data (run SFC, blue). The grey shading comprises the same sensitivity cases. For comparison, the $p_m^{CO_2}$ climatology by Takahashi et al. (2009) is given (violet, with a small additive adjustment to transfer the nominal year 2000 average into a 1997–2009 average according to the rising atmospheric $CO_2$ partial pressure field $p_a^{CO_2}$, Eq. A3).

from scale mismatch between the model pixel ($\approx 4^\circ$ lat. $\times 5^\circ$ long.) and the point measurement, but likely also from the need to compromise between the signals of neighbouring locations in the fit: the statistics of the misfit in the surrounding area (Pacific $45^\circ–90^\circ$ N) reveal very small seasonal biases (well below 1 µatm, Fig. 9). The distribution of residuals deviates somewhat from the Gaussian assumed mathematically in the Bayesian inference by having wider tails, but is symmetric. Example locations STM and – to a lesser degree – BTM show close seasonal fit (Fig. 8), and agreement also in some high-frequency features (see below). A less successful fit is found at Stratus (South Pacific) in an area of low data density (note that the independent data at this location are mainly from after the end of the SOCAT v1.5 data period; thus the agreement relies on the extrapolating effect of the dominating seasonal degrees of freedom (Appendix A2.2))
<figure><img src="image_11.png" /></figure>

**Fig. 8.** Comparison between the estimated $p^{CO_2}$ field (run **SFC**, blue) and measurements (black dots) at test locations. Left: data contained in SOCAT v1.5 (larger dots) at a pixel in the North Pacific, and the pixel around Station M in the North Atlantic. Right: two mooring sites in the North Atlantic (BTM, 32° N, 64° W) and South Pacific (Stratus, 19.70° S, 85.62° W) by Sabine et al. (2010) that are not part of SOCAT v1.5 (smaller dots). The $p^{CO_2}$ field has been picked at the times where the measurements exist (connected by straight lines for clarity); the time axes have been limited to respective years with data. The monthly climatology by Takahashi et al. (2009) (regridded and with added trend parallel to the atmospheric $CO_2$ increase, violet) is also given. See Supplement for comparisons to SOCAT data (Fig. S7.5) and independent data (Fig. S7.6) at further locations.

<figure><img src="image_12.png" /></figure>

**Fig. 9.** Histograms of the residuals of the fit to SOCAT data (i.e. $p^{CO_2}$ estimated by run **SFC** minus data, µatm, horizontal axis). Thin lines show the relative abundance (within bins of 5 µatm) of residuals over all pixels from the North Pacific (north of 45° N) and all time steps with data in winter (January–February–March, blue) or in summer (July–August–September, red). Respective thick lines give a Gaussian of the same mean and standard deviation. See Supplement Fig. S5.2 for global maps of seasonal residuals.

but misses any interannual anomalies). See the Supplement (Sect. S5, Figs. S7.5 and S7.6) for more comparisons and statistics.

In order to validate the high-frequency variations, we compare them to mooring observations independent of the SO-

CAT v1.5 data set (Fig. 10, upper panels). In the extra-tropics (example sites BSO or MOSEAN+WHOTS), many short-term features are reproduced in run **SFC**, though mostly with reduced amplitude. Discrepancies (e.g. beginning of 2007 at BSO) can largely be explained by the deviation of sea surface temperature in our driving field – regridded to the ≈4° × 5° pixels around the stations – from the local temperature (lower panels). Note that this discrepancy does not preclude that the pixel average of $p^{CO_2}$ is actually reproduced more accurately than the local value, to the extent that our SST field represents the pixel’s average temperature; however, this cannot be validated or falsified based on the point data.

In contrast to the extra-tropics, the high-frequency variations in our estimate are unrelated or even opposite to those of the measurements in the tropics (example site TAO170W), even though the SST field performs similarly as at the extra-tropical sites. While in the extra-tropics warmer SST leads to higher $p_m^{CO_2}$ as reproduced by the solubility and chemistry parameterizations, tropical $p_m^{CO_2}$ increases during cold events, as in the upwelling of cold and carbon-rich deep water. Such events would have to be reproduced as high-frequency variations in the ocean-internal sources/sinks $f_{\text{int}}^{\text{DIC}}$, which however cannot be constrained from the SOCAT data set (potentially except at very few locations in periods of high-frequency sampling).

As similar conditions are found at other extra-tropical and tropical sites (Supplement Figs. S7.6, S7.7), we conclude that our results reproduce at least part of the real short-term variations in the extra-tropics dominated by solubility and
<table><caption>**Table 1.** Data sets used as constraints or as driver data in the process parameterizations (Sect. A1).</caption><thead><tr><th>Quantity</th><th>Data set</th><th>Reference</th><th>Pre-treatment/resolution/remarks</th><th>Used for</th></tr></thead><tbody><tr><td colspan="5">Surface-ocean data:</td></tr><tr><td><i>P</i><sub>m</sub><sup>CO2</sup></td><td>SOCATv1.5</td><td>Pfeil et al. (2012), http://www.socat.info/</td><td>Data are used having WOCE-flag = 2 and valid fields for fugacity, temperature, and salinity. Values below 200 µatm or above 600 µatm have been excluded as being local compared to the grid cells. Values have been transferred from fugacity to partial pressure by dividing by 0.996.</td><td>SFC</td></tr><tr><td colspan="5">Atmospheric data:</td></tr><tr><td><i>c</i><sub>CO2</sub></td><td>Various measurement programs</td><td>Conway et al. (1994) and many others</td><td>CO<sub>2</sub> mixing ratios <i>c</i><sub>CO2</sub>. Data points are individual flask pair averages (approximately weekly) or hourly averages, taken at 57 sites globally.</td><td>ATM</td></tr><tr><td colspan="5">Driver data:</td></tr><tr><td><i>T</i></td><td>OAFlux</td><td>Yu and Weller (2007)</td><td>Gridded, daily</td><td>All</td></tr><tr><td><i>ε</i></td><td>OAFlux</td><td>Yu and Weller (2007)</td><td>Gridded, daily</td><td>All</td></tr><tr><td><i>S</i></td><td>WOA 2001</td><td>Conkright et al. (2002)</td><td>Gridded monthly climatology (interpolated, values taken from <i>A</i> data set)</td><td>All</td></tr><tr><td><i>A</i></td><td></td><td>Lee et al. (2006)</td><td>Gridded monthly climatology (interpolated)</td><td>All</td></tr><tr><td><i>γ</i><sup>DIC</sup></td><td></td><td>Egleston et al. (2010)</td><td>Mean spatial pattern</td><td>All</td></tr><tr><td><i>C</i><sub>m</sub><sup>DIC, Ref</sup></td><td>GLODAP</td><td>Key et al. (2004)</td><td>Mean spatial pattern</td><td>All</td></tr><tr><td><i>h</i></td><td>LOCEAN</td><td>de Boyer Montégut et al. (2004)</td><td>Temperature criterion; monthly climatology (interpolated)</td><td>All</td></tr><tr><td><i>u</i></td><td>NCEP reanalysis</td><td>Kalnay et al. (1996)</td><td>Driver for TM3 atmospheric transport model</td><td>All</td></tr><tr><td>Meteo.</td><td>NCEP reanalysis</td><td>Kalnay et al. (1996)</td><td></td><td>All</td></tr><tr><td colspan="5">Comparison data:</td></tr><tr><td><i>P</i><sub>m</sub><sup>CO2</sup></td><td></td><td>Takahashi et al. (2009)</td><td>Gridded monthly climatology</td><td>Sect. 3.5</td></tr><tr><td><i>C</i><sub>m</sub><sup>PO4</sup></td><td>WOA 2005</td><td>Garcia et al. (2006)</td><td>Gridded monthly climatology</td><td>Sect. 3.7</td></tr></tbody></table>

Glossary: GLODAP = Global Ocean Data Analysis Project; LOCEAN = Laboratoire d'océanographie et du climat: expérimentations et approches numériques; Meteo. = Meteorological variables as listed in Heimann and Körner (2003); NCEP = National Centers for Environmental Prediction; OAFlux = Objectively Analysed air-sea Fluxes; SOCAT = Surface Ocean $CO_2$ Atlas; WOCE = World Ocean Circulation Experiment; WOA = World Ocean Atlas; Math symbols see Table 2.

chemistry. In the tropics, more realistic short-term variations could potentially be obtained by parameterization of fast processes (upwelling events) in $f_{int}^{DIC}$.

## 3.5 Comparison to the Takahashi climatology

The monthly mean seasonal cycle calculated from our spatio-temporal $p^{CO_2}$ field has been compared to the widely used observation-based climatology by Takahashi et al. (2009), which differs from our method in the $p^{CO_2}$ data set used (Lamont-Doherty Earth Observatory (LDEO) database, Takahashi et al., 2010) and in the way of spatio-temporal extrapolation of the point data into pixels without data: while the diagnostic scheme extrapolates the internal flux field $f_{int}^{DIC}$ and then infers $p^{CO_2}$ from it via the process representations (Fig. 5), Takahashi et al. (2009) extrapolates the $p^{CO_2}$ field itself, using a 2-dimensional advection-diffusion equation<fnref n="1" />.

In data-dense areas, results should not depend much on the extrapolation method.

In most regions, the seasonal cycle estimated from SOCAT by the diagnostic scheme is similar in phase to the $p^{CO_2}$ climatology by Takahashi et al. (2009) (Fig. 7), though the seasonal amplitude is somewhat larger in many regions. Even better agreement is found when only considering the open ocean: the open-ocean seasonalities in the temperate North Pacific and Atlantic almost fully coincide (Supplement Fig. S7.8). Open-ocean agreement is also confirmed at several of our test locations (Fig. 8).

Differences between our seasonalities and the Takahashi et al. (2009) climatology are somewhat bigger in areas of low data density (e.g. the Southern Ocean), as expected. Comparison to data at typical test locations tends to confirm the slightly larger seasonal amplitudes in our estimate (Fig. 8).

<footnote n="1">¹ As a further methodological difference, the diagnostic scheme uses each data point at its original sampling time; thus no transfer to a nominal year as in Takahashi et al. (2009) is needed and interannual variations in atmospheric $CO_2$ above the ocean are fully taken into account.</footnote>
<table><caption>**Table 2.** Mathematical symbols (also see Table 3). Some quantities are indexed by species or other distinctions explained in the text.</caption><thead><tr><th>Quantity</th><th>Unit<sup>a</sup></th><th>Meaning</th></tr></thead><tbody><tr><td><b>A</b></td><td>ppm (µmol m<sup>-2</sup> s<sup>-1</sup>)<sup>-1</sup></td><td>Atmospheric transport matrix, Eq. (A25)</td></tr><tr><td><b>C</b><sub>m</sub><sup>DIC</sup></td><td>µmol kg<sup>-1</sup></td><td>Mixed-layer DIC concentration</td></tr><tr><td><b>C</b><sub>a</sub><sup>DIC</sup></td><td>µmol kg<sup>-1</sup></td><td>DIC concentration in equilibrium with atmosphere</td></tr><tr><td><b>c</b><sub>obs</sub></td><td>ppm</td><td>Set of observed atmospheric mixing ratios</td></tr><tr><td><b>c</b><sub>mod</sub></td><td>ppm</td><td>Modelled mixing ratios, Eq. (A25)</td></tr><tr><td><b>c</b><sub>0</sub></td><td>ppm</td><td>Initial mixing ratio</td></tr><tr><td><b>f</b><sub>ma</sub></td><td>µmol m<sup>-2</sup> s<sup>-1</sup></td><td>Sea-to-air tracer flux, Eq. (A17)</td></tr><tr><td><b>f</b><sub>hist</sub></td><td>µmol m<sup>-2</sup> s<sup>-1</sup></td><td>History (re-entrainment) flux, Eq. (A20)</td></tr><tr><td><b>f</b><sub>int</sub></td><td>µmol m<sup>-2</sup> s<sup>-1</sup></td><td>Ocean-internal tracer sources/sinks, Eq. (A18)</td></tr><tr><td><b>F</b></td><td>µmol m<sup>-2</sup> s<sup>-1</sup></td><td>Flux model matrix, Eq. (A26)</td></tr><tr><td><b>f</b></td><td>µmol m<sup>-2</sup> s<sup>-1</sup></td><td>Vector of flux discretized at pixels/time steps</td></tr><tr><td><b>h</b></td><td>m</td><td>Mixed layer depth</td></tr><tr><td><b>J</b></td><td>l</td><td>Cost function, Eq. (A27)</td></tr><tr><td><b>J</b><sub>c</sub></td><td>l</td><td>Cost function, atmospheric data part, Eq. (A24)</td></tr><tr><td><b>J</b><sub>p</sub></td><td>l</td><td>Cost function, surface-ocean data part, Eq. (A30)</td></tr><tr><td><b>k</b></td><td>m s<sup>-1</sup></td><td>Piston velocity, Eq. (A2)</td></tr><tr><td><b>L</b></td><td>mol kg<sup>-1</sup> atm<sup>-1</sup></td><td>Solubility (including fugacity correction), Eqs. (A1) and (S1.1)</td></tr><tr><td><b>p</b><sub>m</sub></td><td>µatm</td><td>Partial pressure in mixed layer, Sect. A1</td></tr><tr><td><b>p</b><sub>a</sub></td><td>µatm</td><td>Atmospheric partial pressure, Eq. (A3)</td></tr><tr><td><b>p</b><sub>baro</sub></td><td>atm</td><td>Barometric pressure at ocean surface</td></tr><tr><td><b>p</b><sub>dry</sub></td><td>atm</td><td>Dry-air pressure at ocean surface</td></tr><tr><td><b>p</b><sub>H<sub>2</sub>O</sub></td><td>µatm</td><td>Saturation water vapor pressure above ocean, Eq. (S1.3)</td></tr><tr><td><b>S</b></td><td>%</td><td>Mixed-layer salinity<sup>b</sup></td></tr><tr><td><b>Sc</b></td><td>l</td><td>Schmidt number, Appendix S1.1</td></tr><tr><td><b>T</b></td><td>°C</td><td>Mixed-layer temperature (sea surface temperature, SST)</td></tr><tr><td><b>t</b></td><td>s</td><td>Time coordinate</td></tr><tr><td><b>t</b><sub>i</sub>, <b>t</b><sub>e</sub></td><td>s</td><td>Start time, end time of inversion period</td></tr><tr><td><b>u</b></td><td>m s<sup>-1</sup></td><td>Wind speed at 10 m above surface</td></tr><tr><td><b>x</b></td><td>l</td><td>Vector of adjustable parameters, Eq. (A26)</td></tr><tr><td><b>X</b></td><td>ppm</td><td>Dry-air molar mixing ratio in the atmosphere above the ocean</td></tr><tr><td><b>z</b></td><td>m</td><td>Vertical coordinate</td></tr><tr><td><b>β</b></td><td>l</td><td>Temperature factor of CO<sub>2</sub> partial pressure, Eqs. (A5) and (A6)</td></tr><tr><td><b>Γ</b></td><td>l</td><td>Scaling of piston velocities, Eq. (A2)</td></tr><tr><td><b>γ</b><sup>DIC</sup></td><td>µmol kg<sup>-1</sup></td><td>Buffer factor (response factor), Eq. (A11)</td></tr><tr><td><b>ε</b></td><td>l</td><td>Ice-free fraction of ocean surface (0 = ice-covered, 1 = ice-free)</td></tr></tbody></table>

and Supplement Figs. S7.5 and S7.6). In particular, in the North Pacific where differences between the two methods are biggest, larger seasonalities are confirmed both at the individual pixels (see Fig. 8 for a typical example) and in the larger seasonal bias of about 8.5 µatm in the differences between Takahashi et al. (2009) and SOCAT data points in that area (however, the largest contribution to the difference in the $p^{CO_2}$ average of Fig. 7 comes from the northernmost pixels where data are too sparse for a conclusive validation). Data at Station M (Fig. 8) also confirm larger seasonality in the North Atlantic.

Differences between the two estimates are also bigger towards the coast. Note that part of these coastal differences come from the fact that, in order to cover the entire ocean surface, we extrapolated the Takahashi et al. (2009) clima-

tology using its open-ocean values. For example, this contributes to the differences in the temperate North Atlantic due to the Mediterranean (see Supplement, Fig. S7.5).

Partly, the differences in Fig. 7 can be traced to the fact that the two estimates are based on different $p^{CO_2}$ data sets, though there is considerable overlap between the SOCAT and LDEO databases (this influence has explicitly been tested in the Supplement Sect. S6). Sensitivity to data selection partially arises from interannual variations; in particular, differences in the tropics may be related to the exclusion of El Niño data in Takahashi et al. (2009). Moreover, the synthetic data test (Appendix B) reveals slightly too high seasonality in our results, though most of this mismatch is confined to certain marginal areas (that do not contribute much to the sea-air flux due to low gas exchange).

<footnote><sup>a</sup> ppm = µmol mol<sup>-1</sup>; <sup>b</sup> for the present purposes, salinity in % or on the Practical Salinity Scale are considered interchangeable.</footnote>
<figure><img src="image_13.png" /></figure>

**Fig. 10.** Top: time series of surface-ocean pCO₂ (µatm) at independent mooring sites (not part of SOCATv1.5): Barents Sea Opening (72.5° N, 19.5° E, Arthun et al., 2012), combined MOSEAN and WHOTS moorings (22.77° N, 157.92° W, Sabine et al., 2010), and from the TAO array (2° S, 170° W, Sabine et al., 2010). Comparison between the measurements (black dots) and the pCO₂ field estimated by fitting to the SOCAT data (blue). The estimated field has been picked at the times where measurements exist; one year has been selected for each site. Bottom: sea surface temperature measured simultaneously with pCO₂ (black dots) compared to our SST driving field at this pixel (light grey).

Beyond seasonality, agreement is also found in the im-
plied long-term average global sea-air flux. The 1990–2009
average from run SFC is −1.2 PgC yr⁻¹ (range between
−0.95 PgC yr⁻¹ and −1.45 PgC yr⁻¹ due to sensitivity cases
with lower/higher gas exchange rate, Sect. 3.3). Within un-
certainties, this agrees well with the value of (−1.36 ±
0.6) PgC yr⁻¹ recomputed from Takahashi et al. (2009) by
Wanninkhof et al. (2012, Table 2, “Net flux” and “Coastal
area”) using CCMP winds (Atlas et al., 2011).

In summary, we take the similarity between the seasonaliti-
ties from the two methods based on similar data as an indica-
tion that the estimated pCO2 field in most regions is informed
by the data and not primarily method-dependent.

3.6 Comparison of pCO₂-based and atmospheric
results

Figure 11 compares the SOCAT-based fit from Fig. 6 (run **SFC**) with the fit of the diagnostic scheme to atmospheric CO₂ data (run **ATM**). In the Southern Hemisphere, there is broad agreement in the phasing and amplitude of the seasonal cycle. In particular, the difference between the two runs is much smaller than their differences from the (common) prior, i.e. than the adjustments due to the respective data constraints. Less agreement is found in the tropics and the Northern Hemisphere, in particular in the temperate latitudes.

The ocean fluxes estimated from the atmospheric data
are much less robust (more sensitive to changes in uncer-
tain model details, not shown) than those estimated from the
pCO2 data. As investigated by synthetic data testing (Supple-

ment Sect. S4), this is largely related to the additional degrees
of freedom in adjusting land-to-atmosphere fluxes: signals in
the atmospheric data are partially wrongly attributed to land
or ocean. Consistently, best agreement between **ATM** and
**SFC** occurs in the southern regions away from land influ-
ence, while the large disagreement in the temperate North
Pacific is due to the short but strong sink period during sum-
mer, more typical for terrestrial boreal ecosystems.

The results establish that oceanic and atmospheric data
are consistent in regions away from land influences, but that
estimates based on pCO2 data provide much stronger con-
straints on internal ocean processes, in particular in ocean ar-
eas closer to land. Combining the oceanic and atmospheric
constraints is therefore both possible and beneficial (com-
pare Jacobson et al., 2007). It can be implemented by adding
together the cost function contributions of pCO2 data (from
SFC) and of atmospheric CO2 (from ATM). The estimated
fields in this combined case (not shown) are almost identi-
cal to the SFC run, as expected from the weakness of the
atmospheric constraint on oceanic fluxes. However, as the at-
mospheric data constrain the sum of ocean and land fluxes,
improvements of the estimated land fluxes can be expected.
Of course, the combined run is essentially equivalent to us-
ing the sea-air fluxes of the SFC run as priors in a subsequent
"classical" atmospheric inversion (similar as in the joint in-
version by Jacobson et al., 2007). Details will be discussed
further in the context of interannual variations in the com-
panion paper (Rödenbeck et al., 2013).
<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Monthly mean seasonal cycle of sea-air CO₂ fluxes as Fig. 6. Estimates by fitting the diagnostic scheme to SOCAT data (run **SFC**, blue) are compared to the fit of the scheme to atmospheric CO₂ mixing ratio data (run **ATM**, magenta), as well as to the prior (no data constraint, thin dashed grey)</figcaption></figure>

## 3.7 Plausibility of the ocean-internal carbon sources and sinks

As the diagnostic scheme (in run SFC) is constructed to match the CO₂ partial pressure data, any errors in the parameterizations linking pCO₂ and $f_{\text{int}}^{\text{DIC}}$ (namely carbonate chemistry (including temperature or alkalinity dependence), or mixed-layer carbon budget (including freshwater and history fluxes), see Fig. 1) will be compensated by spurious additional adjustments to the ocean-internal carbon sources and sinks. On the one hand, this means that the estimated pCO₂ (and sea-air flux) fields will be only little affected by imperfections in these parameterizations. Even at pixels not directly constrained by a pCO₂ data value but only indirectly via extrapolation of the $f_{\text{int}}^{\text{DIC}}$ field through its a-priori spatiotemporal correlations (Fig. 5), any such bias will largely cancel out (to the extent that it is the same at the directly constrained and the extrapolated pixels). This leads to the robustness of the pCO₂ estimate shown in Sect. 3.3. Therefore, as

<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Mixed-layer phosphate concentration (monthly climatology, mean subtracted). Seasonality inferred from the SOCAT fit (blue) assumes all ocean-internal carbon sources and sinks are related to phosphate sources and sinks in Redfield proportions. For comparison, the PO₄ climatology from the World Ocean Atlas (WOA, Garcia et al., 2006) is given (green). The grey uncertainty band around the SOCAT-based estimates again gives the range of sensitivity cases as in Fig. 7; note however that calculated PO₄ is also sensitive to several further model elements not included in this range.</figcaption></figure>

long as we are only interested in the pCO₂ (and sea-air flux) fields, $f_{\text{int}}^{\text{DIC}}$ could be just considered a mathematical device in the mapping of pCO₂.

On the other hand, if the ocean-internal fluxes are of interest themselves (or in the envisaged extensions of the scheme using further data streams like oxygen coupled to carbon via $f_{\text{int}}^{\text{DIC}}$), correctness of the parameterizations is needed. On the seasonal timescale considered here, the plausibility of the estimated $f_{\text{int}}^{\text{DIC}}$ field can be tested in the light of mixed-layer phosphate concentrations. As explained in Appendix C, a modelled mixed-layer phosphate field can approximately be calculated from our $f_{\text{int}}^{\text{DIC}}$ estimate, and can be compared to an observation-based monthly climatology (World Ocean Atlas, Garcia et al., 2006). The inferred PO₄ concentration
seasonality is broadly similar in phase and amplitude to the observations (Fig. 12), including characteristic region-to-region differences, though overestimation is seen especially in the temperate latitudes. Given that part of the mismatch arises because the modelled $PO_4$ concentration neglects non-Redfieldian sources/sinks (as justified in Appendix C) and because of errors in the observed $PO_4$ climatology which are likely significant (and amplified by the Redfield ratio), we take this agreement as broad confirmation of the employed parameterizations.

The carbonate chemistry parameterization (including its driving fields, in particular alkalinity) can specifically be checked using DIC data. This is done in the companion paper (Rödenbeck et al., 2013).

## 3.8 The diagnostic mixed-layer scheme in the context of other $p_m^{CO_2}$ mapping methods

The suitability of individual $p_m^{CO_2}$ mapping methods (Sect. 1) depends on the intended application. The primary motivation of the diagnostic mixed-layer scheme presented here was to be able to process different data streams (so far $p^{CO_2}$ and atmospheric $CO_2$ mixing ratio) in a consistent way, and to combine them. The intention was to achieve this in a data-dominated way, with the least complex model possible.

Statistical methods like the reference $p_m^{CO_2}$ climatology by Takahashi et al. (2009) or the interpolation by Jones et al. (2012b) have the advantage of being data-dominated, quite independent of driver data sets and process parameterizations. Regression methods (multi-linear regression, neural networks) use driver data sets to bring in information about modes of variability not resolved by the data, thus addressing the problem of data sparsity; still, they leave flexibility in the way of how $p_m^{CO_2}$ and these driving fields are related. Finally, data assimilation into process models then also prescribes this relation analytically, which brings in process knowledge as a further constraint, with the risk of suppressing real variability not captured in the chosen parameterizations.

The diagnostic scheme is similar to the regression methods, thereby also reproducing some of the short-term variability (Sect. 3.4). It is more rigid, like data assimilation, for some processes that may be considered relatively well known (sea-air gas exchange, carbonate chemistry, and mixed-layer mass conservation). It is more free, on the other hand, for the ocean-internal sources and sinks (it is independent of any assumptions except for smoothness as needed to regularize the mathematics), i.e. it is purely statistical here. The smoothness assumption represents the mechanism of spatio-temporal interpolation<fnref n="2" />.

An advantage of explicit parameterizations, as in the proposed mixed-layer scheme, is that the unknowns (here: internal sources/sinks) are physical quantities which are potentially open for interpretation (compare Sect. 3.7). They also enable us to use measurements of further quantities of the scheme (such as mixed-layer DIC concentration, nutrients, etc.) as data constraints as well. These further data streams can either be used instead of the present ones ($p_m^{CO_2}$ or atmospheric data) to create alternative estimates for comparison, or in conjunction with the present ones to potentially introduce further unknowns and so distinguish more processes (e.g. carbonate chemistry).

All carbon sources/sinks in the ocean and on land are linked to oxygen sinks/sources, most of them in rather well-defined stoichiometric ratios. The diagnostic scheme can therefore be extended to also make use of atmospheric or oceanic measurements of oxygen (Rödenbeck et al., 2013).

Beyond the scheme as presented here, multi-linear regression could be brought in by expressing the unknown ocean-internal fluxes as a linear combination of (further) driving variables and then match the data by adjusting time-independent weights. This would make the estimation more similar to data assimilation, with the benefits and risks mentioned above. Further, instead of linear combinations, parameterizations of the ocean-internal processes involving adjustable parameters could be used if available. This would also open up the way to use even further data streams, such as ocean color from satellite observation.

## 4 Summary and conclusions

We presented a temporally and spatially resolved estimate of the global surface-ocean $CO_2$ partial pressure field and the sea-air $CO_2$ flux, obtained on the basis of $CO_2$ partial pressure data from the SOCAT database. In order to interpolate the point data into a spatio-temporal field and to add short-term variations not constrained from the data, a simple model of surface-ocean biogeochemistry, including representations of sea-air gas exchange, solubility, carbonate chemistry, mixed-layer DIC budget, and seasonal re-entrainment (“history”), was fitted to the observations. The inverse procedure was similar to atmospheric inversion calculations, using spatial and temporal a-priori correlations to extrapolate the data information.

Focusing first on seasonal variations, the following conclusions were drawn:

- The diagnostic scheme can be robustly fit to SOCAT $p^{CO_2}$ data, to estimate the spatio-temporal $p^{CO_2}$ field and sea-air $CO_2$ flux (run SFC).

- In terms of seasonality, the resulting $p^{CO_2}$ field agrees well with the climatology by Takahashi et al. (2009), confirming that the estimates are informed by the data.

<footnote n="2">² This results in different notions of output resolution: Formally, the resolution is daily and $\approx 4^\circ \times 5^\circ$ pixels, but part of this finescale variability is only coming from the driving data, while only the larger scales ($>$ monthly and $>$ 1500 km areas) are actually informed by the SOCAT data.</footnote>
Differences mainly occur in coastal areas and marginal seas, as well as in regions where data density is low.

- The atmospheric $CO_2$ constraint on surface-ocean bio-geochemistry (run **ATM**) is not robust in general, but consistent with the oceanic constraint (run **SFC**) in areas away from land influences. This confirms that it is appropriate and advisable to use the oceanic constraint as oceanic prior in an atmospheric inversion, to improve the inferred land fluxes.

- The diagnostic scheme involves interpretable quantities, as indicated by the comparison of predicted and observed seasonal phosphate variations. The scheme can be extended by parameterizations of ocean-internal sources and sinks, that allow us to use further data streams to constrain individual biogeochemical processes.

- Using synthetic data tests (as in Appendix B), the diagnostic mixed-layer model can also be applied to B assess possible impacts of additional data on constraining ocean biogeochemistry, helping to design future observation strategies.

The results also include interannual variations, which will be considered in a companion paper (Rödenbeck et al., 2013).

The gridded $p^{CO_2}$ and sea-air $CO_2$ flux estimates will be made available for download in digital form from the Jena inversion website, www.bgc-jena.mpg.de. Regular updates are planned.

# Appendix A

## Model documentation

For reference, all mathematical symbols are listed in Tables 2 and 3.

## A1 Ocean process parameterizations

This section details the equations used to express sea-air flux as a function of partial pressure (Sect. A1.1), partial pressure as a function of mixed-layer concentration (Sect. A1.2), and mixed-layer concentration as a function of ocean-interior source/sinks (Sect. A1.3), going down the chain of parameterizations in Fig. 1.

The variables of the following parameterizations are spatio-temporal fields, with a temporal resolution of 1 day and a horizontal resolution of ≈ 4° × 5°, i.e. the equations are applied to every pixel of the atmospheric transport model. Vertically, we consider an oceanic mixed layer assumed to be perfectly homogeneous in temperature and chemical composition, exchanging inorganic carbon with the overlying atmosphere, the underlying ocean interior, neighbouring pixels, and the organic carbon pool.

### A1.1 Solubility and gas exchange

Sea-air flux is expressed in the usual way in terms of a partial pressure difference,

$$f_{\text{ma}}^{\text{CO}_2} = k^{\text{CO}_2} \varrho L^{\text{CO}_2} (p_{\text{m}}^{\text{CO}_2} - p_{\text{a}}^{\text{CO}_2}) \quad \text{(A1)}$$

(atmospheric sign convention: positive = source to atmosphere), where $p_{\text{a}}^{\text{CO}_2}$ is the $CO_2$ partial pressure in the atmosphere, and $p_{\text{m}}^{\text{CO}_2}$ the $CO_2$ partial pressure in the homogeneous mixed layer<fnref n="3" />. The temperature-dependent solubility $L^{\text{CO}_2}$ relates the partial pressure difference to a dissolved $CO_2$ concentration difference. The seawater density $\varrho$ is needed to relate volume-based and mass-based quantities. The diffusive resistance to sea-air gas exchange is described by the piston velocity $k^{\text{CO}_2}$. We use the formulation according to Wanninkhof (1992):

$$k^{\text{CO}_2} = \Gamma u^2 \left( \frac{Sc^{\text{CO}_2}}{Sc^{\text{Ref}}} \right)^{-0.5}, \quad \text{(A2)}$$

calculated from 6-hourly NCEP wind speed $u$ (Kalnay et al., 1996) and Schmidt number $Sc$ given in the Supplement Sect. S1.1. The global scaling factor $\Gamma$ is chosen such that the global mean $CO_2$ piston velocity $k_{\text{Glob}}^{\text{CO}_2}$ matches the value from Table 3 given by Naegler (2009) (average of values by Naegler and Levin (2006), Krakauer et al. (2006), Sweeney et al. (2007), and Müller et al. (2008) inferred from radiocarbon data).

The atmospheric partial pressure

$$p_{\text{a}}^{\text{CO}_2} = X^{\text{CO}_2} p^{\text{dry}} \quad \text{(A3)}$$

inherits variability not only from atmospheric dry-air pressure $p^{\text{dry}}$ (Supplement Sect. S1.1) but also from the atmospheric $CO_2$ molar mixing ratio $X^{\text{CO}_2}$ just above the ocean surface, which has considerable increasing trend, seasonality, and – even over the ocean – synoptic variability. It is taken from the atmospheric $CO_2$ mixing ratio fields provided by the Jena inversion s81_v3.4 (see http://www.bgc-jena.mpg.de/~christian.roedenbeck/download-CO2/). This $X^{\text{CO}_2}$ field has been obtained from a forward run of the atmospheric transport model (TM3, resolution ≈ 4° × 5° × 19 layers), with surface fluxes from a classical atmospheric transport inversion based on atmospheric $CO_2$ measurements. The $X^{\text{CO}_2}$ field constructed that way is compatible with measured $CO_2$ mixing ratios at the sites used, and therefore realistic at least in the rising trend and large-scale seasonality; the sensitivity of the results presented here to the $X^{\text{CO}_2}$ field is small (see Supplement Sect. S2).<fnref n="4" />

<footnote n="3">³ More precisely, it is the $CO_2$ partial pressure in air in equilibrium with sea water of the present chemical composition.</footnote>

<footnote n="4">⁴ Of course, the use of inversion-based $X^{\text{CO}_2}$ already creates some dependence of our prior from atmospheric data. In run **ATM**, this formally violates assumptions of Bayesian inference, which however could only affect calculations of a-posteriori covariances. Even there, we would not expect problems due to the small sensitivity on $X^{\text{CO}_2}$.</footnote>
<table><caption>**Table 3.** Physical constants, and geophysical quantities assumed constant.</caption><thead><tr><th>Quantity</th><th>Value</th><th>Meaning</th></tr></thead><tbody><tr><td>c<sub>p</sub></td><td>3993 J kg<sup>-1</sup> K<sup>-1</sup></td><td>Specific heat capacity of sea water (value at 20 °C)</td></tr><tr><td>k<sup>CO<sub>2</sub></sup><sub>Glob</sub></td><td>16.5 cm h<sup>-1</sup></td><td>(= 4.58 · 10<sup>-5</sup> m s<sup>-1</sup>) Global mean piston velocity of CO<sub>2</sub> (Naegler, 2009)</td></tr><tr><td>R<sub>gas</sub></td><td>8.3144 J mol<sup>-1</sup> K<sup>-1</sup></td><td>Gas constant</td></tr><tr><td>T<sub>0</sub></td><td>273.15 K</td><td>Absolute temperature at 0 °C</td></tr><tr><td>V<sub>ideal</sub></td><td>0.0224136 m<sup>3</sup> mol<sup>-1</sup></td><td>Molar volume of an ideal gas</td></tr><tr><td>ρ</td><td>1025 kg m<sup>-3</sup></td><td>Density of sea water</td></tr></tbody></table>

<table><caption>**Table 4.** Adjustable degrees of freedom.</caption><thead><tr><th>Quantity</th><th></th><th>A-priori value</th><th>Spatial resolution</th><th>Temporal resolution</th></tr></thead><tbody><tr><td colspan="5">Carbon cycle components:</td></tr><tr><td>f<sup>DIC</sup><sub>int</sub></td><td>Ocean-internal DIC sources/sinks</td><td>Zero</td><td>corr. pix.</td><td>lt + seas</td></tr><tr><td>f<sup>CO<sub>2</sub></sup><sub>nee</sub></td><td>Land NEE (run <b>ATM</b> only)</td><td></td><td>corr. pix.</td><td>lt, seas, var</td></tr><tr><td colspan="5">Technical degrees of freedom:</td></tr><tr><td>ΔC<sup>DIC, ini</sup><sub>m</sub></td><td>Initial condition of budget equation</td><td>Zero</td><td>corr. pix.</td><td>–</td></tr><tr><td>f<sup>CO<sub>2</sub></sup><sub>ini</sub></td><td>Flux pulse for atmospheric initial condition (run <b>ATM</b> only)</td><td>Zero</td><td>3 lat. bands</td><td>–</td></tr></tbody></table>

Glossary: lt, – = constant (long-term); seas = seasonal; var = interannual and short-term variability; corr. pix. = correlated pixels; math symbols see Table 2.

## A1.2 Carbonate chemistry

The CO₂ partial pressure in the mixed layer is a function of mixed-layer DIC concentration $C_m^{DIC}$, alkalinity A, temperature T, and salinity S:

$$p_m^{CO_2} = p_m^{CO_2} (C_m^{DIC}, A, T, S). \qquad (A4)$$

In the diagnostic scheme, the relevant dependence is that on $C_m^{DIC}$, while the influences of A, T, and S are calculated a-priori from driver data (Table 1). While the exact function is a 5th order polynomial (Zeebe and Wolf-Gladrow, 2001), we use existing approximations, as described in the following.

At first, we single out the temperature dependence into a temperature-dependent factor and $p_m^{CO_2}$ at a reference temperature $T^{Ref}$:

$$p_m^{CO_2} = \beta(T, T^{Ref}) \cdot p_m^{CO_2}(C_m^{DIC}, A, T^{Ref}, S). \qquad (A5)$$

For iso-chemical sea water (i.e. constant $C_m^{DIC}$, A, and S), Takahashi et al. (2009) experimentally determined the ratio of $p_m^{CO_2}$ of samples at different temperatures to be

$$\beta(T_1, T_2) = \exp \left\{ \begin{array}{l} 0.0433 \text{ K}^{-1} (T_1 - T_2) \\ -4.35 \times 10^{-5} \text{ K}^{-2} (T_1^2 - T_2^2) \end{array} \right\}. \qquad (A6)$$

Even though this relation has been fitted to data in the temperature range of about 2 to 25 °C, it is used here for all temperatures.

The remaining dependences of $p_m^{CO_2}$ are non-linear, but monotonic and can be linearized to good approximation (Sarmiento and Gruber, 2006). Linearization of the dependence on $C_m^{DIC}$ is important here in order to be able to use the fast minimization algorithm of Rödenbeck (2005) (Sect. A2.1). We consider deviations ($\Delta$) of the driver variables from temporally constant reference values (superscript “Ref”, see Sect. A1.4 below):

$$C_m^{DIC} = \Delta C_m^{DIC} + C_m^{DIC, Ref} \qquad (A7)$$

$$A = \Delta A + A^{\text{Ref}} \qquad (A8)$$

$$S = \Delta S + S^{\text{Ref}} \qquad (A9)$$

Linearization of Eq. (A5) around these references gives

$$p_m^{CO_2} = \beta(T, T^{Ref}) \times \left[ p_m^{CO_2, Ref} + \frac{\partial p_m^{CO_2}}{\partial C_m^{DIC}} \Delta C_m^{DIC} + \frac{\partial p_m^{CO_2}}{\partial A} \Delta A + \frac{\partial p_m^{CO_2}}{\partial S} \Delta S \right] \qquad (A10)$$

The DIC sensitivity has been calculated from the response factors $\gamma^{DIC}$ by Egleston et al. (2010):

$$\frac{\partial p_m^{CO_2}}{\partial C_m^{DIC}} = \frac{p_m^{CO_2, Ref}}{\gamma^{DIC}}. \qquad (A11)$$
The alkalinity sensitivity has been calculated using the approximation

$$ \frac{\partial p_m^{CO_2}}{\partial A} = \frac{\partial p_m^{CO_2, \text{Ref}}}{A^{\text{Ref}}} \times \frac{-(A^{\text{Ref}})^2}{(2C_m^{\text{DIC, Ref}} - A^{\text{Ref}})(A^{\text{Ref}} - C_m^{\text{DIC, Ref}})} \quad (A12) $$

(Eq. (8.3.17) of Sarmiento and Gruber, 2006, note missing minus sign). The salinity sensitivity (direct effect on chemical equilibrium only) has been set to

$$ \frac{p_m^{CO_2, \text{Ref}}}{\partial S} = \frac{p_m^{CO_2, \text{Ref}}}{S^{\text{Ref}}} \quad (A13) $$

(Eq. (8.3.5) of Sarmiento and Gruber, 2006). Any seasonal and interannual variations in the proportionalities $\partial p_m^{CO_2}/\partial C_m^{DIC}$, etc., are neglected. The variations in alkalinity, A, and salinity, S, are fixed according to input data sets (Table 1). Their effect on $p_m^{CO_2}$ is illustrated in the Supplement Fig. S7.2. Note that the effect of the alkalinity variations, which are strongly related to freshwater fluxes, is mainly offset by the freshwater effects on DIC considered below (Eq. A19).

Equation (A10) provides the desired linear link between changes in $p_m^{CO_2}$ and $C_m^{DIC}$ according to carbonate chemistry. For simplification in the further use, we combine the gas exchange and chemistry formulas: we first define an apparent DIC concentration in equilibrium with the atmosphere,

$$ \Delta C_a^{\text{DIC}} = \left( \frac{\partial p_m^{\text{CO}_2}}{\partial C_m^{\text{DIC}}} \right)^{-1} \left[ \frac{p_a^{\text{CO}_2}}{\beta(T, T^{\text{Ref}})} - p_m^{\text{CO}_2, \text{Ref}} \right. \quad (A14) \\ \left. - \frac{\partial p_m^{\text{CO}_2}}{\partial A} \Delta A - \frac{\partial p_m^{\text{CO}_2}}{\partial S} \Delta S \right] $$

Using this definition, Eq. (A10) becomes

$$ p_m^{\text{CO}_2} = p_a^{\text{CO}_2} + \beta(T, T^{\text{Ref}}) \frac{\partial p_m^{\text{CO}_2}}{\partial C_m^{\text{DIC}}} (\Delta C_m^{\text{DIC}} - \Delta C_a) \quad (A15) $$

Now also defining an “apparent piston velocity” of DIC as

$$ k^{\text{DIC}} = k^{\text{CO}_2} \cdot L^{\text{CO}_2} \beta(T, T^{\text{Ref}}) \frac{\partial p_m^{\text{CO}_2}}{\partial C_m^{\text{DIC}}} \quad (A16) $$

Eqs. (A1) and (A15) can be combined and simplified into

$$ f_{\text{ma}}^{\text{CO}_2} = k^{\text{DIC}} \varrho (\Delta C_m^{\text{DIC}} - \Delta C_a), \quad (A17) $$

expressing the sea-air $CO_2$ flux as a linear function of the surface DIC concentration.

## A1.3 Mixed-layer DIC budget

To finally establish the link between the mixed-layer DIC concentration and the ocean-interior sources and sinks, we

consider that changes in the spatio-temporal field of dissolved carbon in the mixed layer need to be balanced by the sum of fluxes:

$$ \frac{d}{dt} C_m^{\text{DIC}} = \frac{1}{h\varrho} \left( -f_{\text{ma}}^{\text{CO}_2} + f_{\text{hist}}^{\text{DIC}} + f_{\text{frw}}^{\text{DIC}} + f_{\text{int}}^{\text{DIC}} \right). \quad (A18) $$

Fluxes (carbon exchange per unit horizontal area and unit time) comprise loss through sea-air exchange $f_{\text{ma}}^{\text{CO}_2}$, seasonal re-entrainment $f_{\text{hist}}^{\text{DIC}}$ and freshwater effects $f_{\text{frw}}^{\text{DIC}}$ explained below, and fluxes $f_{\text{int}}^{\text{DIC}}$ due to all other ocean-internal processes (biological conversion, vertical advection/diffusion, or horizontal advection/diffusion). This formulation assumes that fluxes are immediately diluted vertically throughout the mixed layer. Mixed-layer depth h is prescribed as a climatology by de Boyer Montégut et al. (2004) (downloaded from http://www.locean-ipsl.upmc.fr/~cdblod/mld.html using the MLD climatology mldT02-sk.nc according to a temperature criterion; chosen because the MLD climatology based on a density criterion has incomplete coverage, and there is no MLD product directly based on data that also represents interannual variations). The density $\varrho$ of surface water (assumed constant, Table 3) is needed as $C_m^{\text{DIC}}$ is mass based (given in µmol kg$^{-1}$), while fluxes are volume/area based.

Two specific ocean-internal processes are not considered part of $f_{\text{int}}^{\text{DIC}}$ but are parameterized explicitly<fnref n="5" /> in the budget (Eq. A18):

- The freshwater effect $f_{\text{frw}}^{\text{DIC}}$ describes the dilution or enhancement of the mixed-layer DIC concentration by precipitation, evaporation, river freshwater input, or sea-ice formation/melting. Analogously to the often used salinity normalization (Sarmiento and Gruber, 2006), we assume that freshwater does not contain either carbon or salt (and that all salinity variations are related to freshwater fluxes), and calculate the freshwater effect on DIC from salinity changes:

$$ f_{\text{frw}}^{\text{DIC}} = h\varrho \frac{C_m^{\text{DIC,Ref}}}{S^{\text{Ref}}} \frac{dS}{dt}. \quad (A19) $$

This will tend to overestimate the effect to the extent that freshwater input does carry carbon.

- The “history” flux $f_{\text{hist}}^{\text{DIC}}$ is part of the carbon exchange between the mixed layer and the underlying ocean. Consider the hypothetical case that no sources or sinks are acting below the mixed layer. Then, as soon as the mixed layer is deepening, it re-entrains water that was left behind previously during mixed-layer shoaling. As mixed-layer concentrations will generally have changed

<footnote n="5">⁵ The explicit treatment of these processes is not strictly needed if the only purpose of the scheme was mapping of $p^{CO_2}$, because they would otherwise be absorbed into $f_{\text{int}}^{\text{DIC}}$ during the fit to the data. However, splitting them off allows easier interpretation of $f_{\text{int}}^{\text{DIC}}$ (Sect. 3.7), and also prepares for the envisaged link to oxygen.</footnote>
meanwhile, this re-entrainment is equivalent to a carbon flux which is proportional to the difference between the current concentration $C_m^{\text{DIC}}(t)$ and that at time $t_{\text{prev}}$ when the mixed layer was as deep for the last time. We call this the “history” flux. It can be written as

$$f_{\text{hist}}^{\text{DIC}}(t) = \varrho \cdot (C_{\text{m}}^{\text{DIC}}(t_{\text{prev}}) - C_{\text{m}}^{\text{DIC}}(t)) \cdot \vartheta \left(\frac{dh}{dt}\right) \quad (\text{A20})$$

The Heaviside step function $\vartheta$ ensures that $f_{\text{hist}}$ only acts during mixed-layer deepening ($dh/dt > 0$). The history flux is non-zero in the long-term mean: The deeper waters that get disconnected from the shoaling mixed layer will preserve the high DIC concentrations of the winter season. A budget equation that does not include the re-entrainment of this amount of DIC during mixed-layer deepening would imply a systematic downward “leak” of DIC. The history flux is thus essential to balance the long-term mean internal sources/sinks $f_{\text{int}}^{\text{DIC}}$ with the mean sea-air exchange $f_{\text{ma}}^{\text{CO}_2}$ (plus the mean DIC accumulation in the mixed layer itself). This has been verified numerically, and Sect. S3 proves this balance analytically.

Any changes of the DIC concentration below the mixed layer will lead to additional entrainment fluxes, forming an important part of $f_{\text{int}}^{\text{DIC}}$.

Substituting Eq. (A17) into Eq. (A18) gives

$$\frac{d}{dt} \Delta C_{\text{m}}^{\text{DIC}} = -\frac{k^{\text{DIC}}}{h} \Delta C_{\text{m}}^{\text{DIC}} + \frac{k^{\text{DIC}}}{h} \Delta C_{\text{a}}^{\text{DIC}} + \frac{1}{h\varrho} f_{\text{hist}}^{\text{DIC}} \{\Delta C_{\text{m}}^{\text{DIC}}\} + \frac{1}{h\varrho} f_{\text{frw}}^{\text{DIC}} + \frac{1}{h\varrho} f_{\text{int}}^{\text{DIC}}, \quad (\text{A21})$$

where we further exploited that the temporally constant reference $C_m^{\text{DIC, Ref}}$ disappears from the left-hand side. The reference also cancels out in the history flux (Eq. A20), which is thus a linear functional of $\Delta C_m^{\text{DIC}}(t)$. Equation (A21) thus represents a 1st order differential equation in $\Delta C_m^{\text{DIC}}$ that can be solved (numerically) for given “forcing” $f_{\text{int}}^{\text{DIC}}$. At the beginning $t_i$ of the time period, we use the initial condition

$$\Delta C_{\text{m}}^{\text{DIC}}(t_i) = \Delta C_{\text{a}}(t_i) + \Delta C_{\text{m}}^{\text{DIC, ini}}. \quad (\text{A22})$$

Up to a constant offset (determined by $C_m^{\text{DIC, Ref}}$) and an initial transient (determined by $\Delta C_m^{\text{DIC, ini}}$, see Supplement Sect. S1.2), Eq. (A21) thus provides a linear dependence of $C_m^{\text{DIC}}$ on $f_{\text{int}}^{\text{DIC}}$ for every pixel.

### A1.4 Reference values

For the reference values $C_m^{\text{DIC, Ref}}$, $A^{\text{Ref}}$, $T^{\text{Ref}}$, and $S^{\text{Ref}}$, we use the long-term averages of the respective driver data sets (Table 1) for each pixel. The partial pressure reference value $p_m^{\text{CO}_2, \text{Ref}}$ in Eq. (A10) formally follows from the functional relation

$$p_m^{\text{CO}_2, \text{Ref}} = p_m^{\text{CO}_2} (C_m^{\text{DIC, Ref}}, A^{\text{Ref}}, T^{\text{Ref}}, S^{\text{Ref}}). \quad (\text{A23})$$

Solving this equation would be numerically involved but fortunately is not needed in practice: Firstly, $C_m^{\text{DIC, Ref}}$ only determines the absolute levels of DIC (Eq. A7) not considered here, but does not directly affect our target quantities $p_m^{\text{CO}_2}$ or sea-air fluxes (as it cancels out both from Eqs. (A15), (A17), and (A21); we nevertheless use a data-based field because a minor influence of $C_m^{\text{DIC, Ref}}$ arises through Eq. (A12) (alkalinity sensitivity) and Eq. (A19) (freshwater effect)). Secondly, the partial pressure reference value $p_m^{\text{CO}_2, \text{Ref}}$ has no effect either: Though changes in $p_m^{\text{CO}_2, \text{Ref}}$ shift the long-term mean of $\Delta C_a$ through Eq. (A14), the same shift then happens to $\Delta C_m^{\text{DIC}}$ through Eq. (A21) with Eq. (A22), such that it again cancels out from Eqs. (A15) and (A17). The choice of $p_m^{\text{CO}_2, \text{Ref}}$ is thus arbitrary (we use $p_a^{\text{CO}_2}(t_i)$ as for numerical reasons a value on the order of the actual $p_m^{\text{CO}_2}$ is desirable).

### A2 Inversion

The estimation procedure of this paper is derived from the linear Bayesian atmospheric inversion (Newsam and Enting, 1988). The implementation is based specifically on the atmospheric transport inversion (“Jena CO₂ inversion”) described in Rödenbeck (2005) and reviewed in Sect. A2.1. This atmospheric transport inversion is extended by representing sea-air fluxes by the ocean parameterizations (Sect. A2.2) and by using the $p^{\text{CO}_2}$ constraint instead of or in addition to the atmospheric constraint (Sect. A2.3).

#### A2.1 The classical atmospheric transport inversion

The atmospheric inversion is based on a set of observed mixing ratios $c_{\text{obs}}$ (time series at observation sites, see Table 1). The inversion calculation seeks fluxes $f$ that lead to the best match between $c_{\text{obs}}$ and corresponding modelled mixing ratios $c_{\text{mod}}(f)$, in the sense that the value of the cost function

$$J_c = \frac{1}{2} (c_{\text{obs}} - c_{\text{mod}})^{\text{T}} Q_c^{-1} (c_{\text{obs}} - c_{\text{mod}}) \quad (\text{A24})$$

is minimal (superscript T means vector transpose). The (diagonal) matrix $Q_c$ introduces a weighting among the concentration values, involving assumed measurement uncertainty, location-dependent model uncertainty, and a data density weighting, as described in Rödenbeck (2005). The modelled mixing ratios, taken at the same time and location as the observations, are simulated by a numerical transport model (here, the TM3 model, Heimann and Körner, 2003), which can formally be written as

$$c_{\text{mod}} = \mathbf{A}f + c_0 \quad (\text{A25})$$

with the transport matrix $\mathbf{A}$ and the initial concentration $c_0$.

Since the atmospheric data at the discrete set of sites cannot fully constrain the flux field at all pixels and time steps, Bayesian a-priori constraints are introduced to regularize the
estimation. They are implemented here by writing the fluxes as a function of a set of parameters x:

$$f(\mathbf{x}) = f_{\text{pri}} + \mathbf{F}\mathbf{x}. \qquad (\text{A26})$$

In this linear flux model (Rödenbeck, 2005), each parameter in x acts as a multiplier for one of the columns of the matrix F to be explained below. A-priori, the parameters x are assumed to have zero mean and unit variance and to be uncorrelated, expressed by adding the Bayesian cost function contribution

$$J = J_c + \frac{1}{2}\mathbf{x}^T\mathbf{x}. \qquad (\text{A27})$$

This is equivalent to assuming fluxes with a-priori mean $f_{\text{pri}}$ and a-priori covariance matrix $Q_f$, $\text{pri} = \mathbf{F}\mathbf{F}^T$.

The specification of the flux model elements $f_{\text{pri}}$ and F is detailed in Rödenbeck (2005). The flux is first written as a sum of contributions from several components:

$$f = f_{\text{ma}} + f_{\text{nee}} + f_{\text{fos}} + f_{\text{ini}} \qquad (\text{A28})$$

$(f_{\text{ma}} = \text{sea-air flux}, f_{\text{nee}} = \text{terrestrial net ecosystem exchange (NEE)}, f_{\text{fos}} = \text{fossil fuel burning emissions}, f_{\text{ini}} = \text{flux pulse at beginning of inversion period to adjust initial atmospheric mixing ratio field). Each component } i \text{ is represented by a structure as Eq. (A26), i.e.}$

$$f(\mathbf{x}) = \sum_{i} (f_{\text{pri},i} + \mathbf{F}_{i}\mathbf{x}_{i}), \qquad (\text{A29})$$

where the vectors $\mathbf{x}_i$ denote subsets of parameters in x, and the matrices $\mathbf{F}_i$ denote the corresponding groups of columns of F. The columns of $\mathbf{F}_i$ represent spatio-temporal base functions out of which the spatially and temporally varying adjustments to the flux field $f_i$ are composed. In space, the base functions represent flux elements localized at each of the grid cells of the transport model. These elements overlap each other by exponential tails, acting to smooth the flux field. In time, the base functions represent Fourier modes. Smoothing is achieved by down-weighting the higher-frequency modes. In addition to the a-priori correlations implemented that way, all columns of $\mathbf{F}_i$ are proportional to a spatio-temporal weighting allowing flux adjustments in areas of activity of the component (e.g. over the ocean for $f_{\text{ma}}$) while suppressing flux adjustments elsewhere.

The detailed flux model settings for the components $f_i$ of the surface-to-atmosphere CO₂ flux in the “Jena inversion” are

- $f_{\text{ma}}$: sea-air flux, being of central interest here. In the pure transport inversion as in Rödenbeck (2005), this component is estimated directly by decomposition into flux elements as just described.

- $f_{\text{nee,lt}} + f_{\text{nee,seas}} + f_{\text{nee, var}}$: terrestrial net ecosystem exchange (NEE), further split into long-term flux

(superscript “lt”), mean seasonality (“seas”), and non-seasonal variations (“var”, interannual and high-frequency variations). The specification of the priors, a-priori sigmas, and correlation structure is taken from the Jena inversion v3.2 (as documented in Rödenbeck (2005) except for a-priori uncertainties tightened by the factor $\sqrt{8}$, and length scales of a-priori spatial correlations increased by a factor 3 in longitude direction and by 1.5 in latitude direction).

- $f_{\text{fos}}$: fossil fuel burning. Only a fixed (prior) term, as in the Jena inversion v3.2 (yearly emissions from the EDGAR v4.0 database, update of Olivier and Berdowski, 2001).

- $f_{\text{ini}}$: purely technical component. Flux pulse at beginning of inversion period to adjust initial atmospheric mixing ratio field (see Rödenbeck, 2005).

The numerical minimization of the cost function is done by Conjugate Gradient descent with re-orthonormalization (Rödenbeck, 2005).

## A2.2 Extension: using ocean parameterization

The inversion set-up used in this study is identical to the Jena inversion system as of Sect. A2.1 in most respects, except that the sea-air flux $f_{\text{ma}}$ is not estimated directly but implemented as a function of the ocean-internal flux $f_{\text{int}}$ using the parameterizations of Sect. A1. Then, instead of $f_{\text{ma}}$, the ocean-internal flux $f_{\text{int}}$ is adjusted to match the data constraint. The detailed flux-model specifications for $f_{\text{int}}$ are similar to those for the sea-air flux in the atmospheric inversion:

- The a-priori state has been defined by $f_{\text{int, pri}}^{\text{DIC}} = 0$ (no internal sources or sinks). The “error” of this prior model $f_{\text{int, pri}}^{\text{DIC}} = 0$ is then identical with the process fluxes themselves. Therefore, the a-priori uncertainty around this prior should reflect the expected amplitude of variations. In the standard set-up, we make no assumptions about spatial structure of the amplitude, and take a constant per-ocean-area uncertainty over all the ocean. No adjustments are allowed in ice-covered regions (see Supplement S1.3). The a-priori uncertainties are scaled such that the implied a-priori uncertainty of the globally integrated sea-air flux is $1200 \text{ Tmol yr}^{-1}$ (with respect to 3-monthly anomalies).

- Bayesian a-priori spatial and temporal correlations have been implemented that dampen variability in $f_{\text{int}}^{\text{DIC}}$ on scales smaller than around 640 km (longitude), 320 km (latitude), and about 2 weeks (time), Fig. A1. The chosen e-folding lengths are one-third of the ocean flux correlations in the standard Jena Inversion v3.4 (update of
<figure><img src="image_16.png" /></figure>

**Fig. A1.** Top: map of a-priori correlation coefficients with respect to an example pixel at (70° E, 49° S), illustrating how inverse adjustments to $f_{\text{int}}^{\text{DIC}}$ de-correlate with distance. Bottom: a-priori correlation coefficients with respect to an example time step in the middle of the inversion period; in addition to the de-correlation away from the example time step, correlations rise again each year due to the dominance of the seasonal cycle.

Rödenbeck, 2005)<fnref n="6" />, and similar to those of pCO₂ found from correlation analysis along cruise tracks by Jones et al. (2012a). The correlations are seen as a mathematical device to stabilize the estimation by suppressing noise which results, for example, from unequal sampling (they represent the main mechanism by which the information from the discrete atmospheric or oceanic data points is extrapolated into space and time, see Fig. 5 for illustration). The correlation scales need to be large enough to accomplish this extrapolation, but small enough not to unduly dampen actual signals. Correlations shorter than in the standard atmospheric inversion have been chosen to reflect the higher spatial detail in the pCO₂ data. For example, this is important in the North Atlantic where shorter correlations are needed to prevent a large influence of coastal data on

the regional average (test not shown). Shorter correlations also prevent implausible leakage of interannual variability of the tropical Pacific northward and southward (see Rödenbeck et al., 2013). On the other hand, longer correlations may have been beneficial in the data-sparse Southern Ocean as they lead to slightly better performance in synthetic-data tests (analogously to Appendix B, test not shown). In any case, simple distance-based de-correlation necessarily represents a crude approximation to the unknown real correlation structure of the ocean-internal sources/sinks determined by complex processes. To ensure that this approximation does not critically determine the result, runs with different de-correlation lengths are included in the set of sensitivity cases in our robustness test (Sect. 3.3).

In the time domain, we additionally dampen all non-seasonal variability by a factor 16 with respect to the seasonal variability (this is accomplished by implementing the time dependence of $f_{\text{int}}^{\text{DIC}}$ as a Fourier series and reducing the a-priori uncertainty of the non-seasonal

<footnote n="6">⁶ As an additional minor technical change, the basic spatial grid of the adjustable degrees of freedom is every pixel of the TM3 transport model grid, rather than aggregates of 4 neighbouring pixels, reflecting the higher spatial resolution of the data.</footnote>
modes, see Rödenbeck, 2005). Such tighter a-priori uncertainties for the interannual variations are consistent with the expectation that their amplitudes are smaller than those of the seasonal variations. Moreover, they dampen spurious interannual variations that may otherwise arise from the very different coverage of the SOCAT data in the individual years, and extrapolate the mean seasonal cycle of $f_{int}^{DIC}$ throughout the inversion period. As the chosen dampening factor is largely arbitrary, it has also been varied in sensitivity tests.

A further degree of freedom is $\Delta C_m^{DIC, ini}$, needed only for technical reasons: Through adjusting $\Delta C_m^{DIC, ini}$, the initial conditions of the carbon budget equation can become consistent with the data (see Supplement Sect. S1.2). Adjustments to $\Delta C_m^{DIC, ini}$ are spatially resolved with the same spatial correlations as $f_{int}$, but constant in time; the a-priori value is zero.

In run **ATM**, there are further degrees of freedom not needed in run **SFC**: as in the atmospheric inversion, land fluxes are adjustable as specified in Sect. A2.1, as well as the adjustable atmospheric initial condition.

Table 4 summarizes the variables that are adjusted in the cost function minimization.

### A2.3 Extension: the $p^{CO_2}$ constraint

Analogously to the observed atmospheric mixing ratios $c_{obs}$, all selected partial pressure values from the SOCAT database (Table 1) are collected into a vector $p_{obs}$. A corresponding vector $p_{mod}$ of modelled partial pressures is formed by sampling the gridded $p_{m}^{CO_2}$ field of the model at the locations and times of the measurements (i.e. the modelled value is that found in the gridcell/timestep which encloses the location/time of the corresponding observation). If several SOCAT data points fall into the same gridcell and timestep, they are averaged together and only form a single element in $p_{obs}$ and $p_{mod}$.

The vector $p_{mod}$ is a function of the adjustable variable $f_{int}$, as is $p_{mod}$ (except that $p_{mod}$ is also a function of further adjustable variables, $f_{nee}$ and $f_{ini}$, which are irrelevant now). We can define an analogous cost function contribution:

$$J_p = \frac{1}{2} (\mathbf{p}_{\text{obs}} - \mathbf{p}_{\text{mod}})^{\text{T}} \mathbf{Q}_p^{-1} (\mathbf{p}_{\text{obs}} - \mathbf{p}_{\text{mod}}). \quad (\text{A30})$$

The covariance matrix $Q_p$ is chosen diagonal as well; the uncertainty for every individual pixel with two or more SOCAT data points is set to 10 µatm (for pixels with a single data point to $\sqrt{2} \times 10$ µatm). The fit to the SOCAT data is then done in the same way as to the atmospheric data, just with $J_c$ replaced by $J_p$ (those degrees of freedom that do not depend on $p_{m}^{CO_2}$ – land fluxes and atmospheric initial condition – are unconstrained in this inversion and just remain at their a-priori values).

For the case combining atmospheric and oceanic constraints (Sect. 3.6), both cost function contributions are added ($J = J_c + J_p + \frac{1}{2}x^T x$).

## Appendix B

### Test of retrieval capability

As a prerequisite to interpreting the results of fitting the scheme to surface-ocean $p^{CO_2}$ data, we need to test whether (1) the scheme contains sufficient degrees of freedom to reproduce the spatial and temporal variability of the $p^{CO_2}$ field, and (2) the information available in the SOCAT data set is sufficient to constrain it. This test is done by a synthetic run: (1) We create a spatio-temporal $p^{CO_2}$ field to be used as synthetic “truth” on the basis of a biogeochemical process model simulation (NEMOv2.3 with PlankTOM5, Buitenhuis et al., 2010), mapped to the grid of our scheme. The seasonal cycle of this field is similar to the observations at most locations except for the high latitudes, which is why we subtract a mean seasonal cycle and replace it by that of Takahashi et al. (2009). This combined field then contains variations similar to reality on all relevant timescales (though variations on the fastest (weekly to daily) timescales are smaller than observed). (2) This synthetic $p^{CO_2}$ field is sampled at the locations and times of the SOCAT measurements, in the same way as the modelled $p^{CO_2}$ field (Sect. A2.3). This gives a set of pseudo data representing the same amount of information (in terms of data density available to the scheme) as the actual data (note that SOCAT also contains along-track variability on scales much shorter than the size of our grid cells, which the scheme cannot make use of in any case). (3) We then fit the diagnostic scheme to the pseudo data, and compare the result with the synthetic truth as the known correct answer.

In terms of mean seasonality considered here, the retrieval (blue) fits the “known truth” (violet) closely (Fig. B1). Most region-to-region differences are reproduced. The largest differences are found in the North Pacific, and in parts of the temperate Southern Hemisphere. If $p^{CO_2}$ is averaged over open ocean only (not shown), Northern Hemisphere and tropical mismatches further reduce, in particular with the essentially perfect match in the temperate North Atlantic.

To be able to judge the match, we set it into perspective with two further cases representing situations with less or with more available information, respectively:

- The black line shows the Bayesian prior, which does not yet use any SOCAT or atmospheric data (Sect. 2.3). The $p^{CO_2}$ seasonality thus just represents a-priori knowledge on the response to temperature-induced changes in the chemistry (and to lesser extend in solubility, Sect. A2.2). It is opposite in phase compared to the synthetic truth in high latitudes. The pseudo data are able to correct this large difference almost completely.
<figure><img src="image_17.png" /><figcaption>Fig. B1. Testing the capacity to retrieve the $p^{CO_2}$ field from the SOCAT data set: fit of the diagnostic scheme to the $p^{CO_2}$ field simulated by the PlankTOM5 biogeochemical model (Buitenhuis et al., 2010) (combined with seasonality from Takahashi et al. (2009)) sub-sampled at locations/times of the SOCAT measurements (blue), compared to the simulated $p^{CO_2}$ field itself (“known truth”, violet). For further comparison, a fit of the scheme to the simulation at every grid cell and time step (“complete information”, green) and the Bayesian prior (“no data information”, thin dashed grey, partially off-scale) is given. The $p^{CO_2}$ fields have been averaged over the regions as in Fig. 7.</figcaption></figure>

- We also did the retrieval based on maximum information about the $p^{CO_2}$ field available, by constraining the scheme at every grid point and time step, not only where actual SOCAT data exist. The result (green) fits the known truth almost perfectly, showing that there are sufficient degrees of freedom in the scheme to reproduce the mean seasonality. However, as the performance of the SOCAT-sampled pseudo data is only slightly worse, we conclude that the SOCAT data density is essentially sufficient to constrain the seasonality of surface-ocean biogeochemistry on the scale of the considered regions (note that strong modes of variability existing in real-

ity but not contained in our synthetic “truth” potentially deteriorate the performance of our fit to real data).

The test confirms that the differences between our results and Takahashi et al. (2009) in the northernmost part of the North Pacific are due to the lack of data (see Supplement Fig. S7.4) in conjunction with spatial variability: In Takahashi et al. (2009), seasonal amplitudes are large in a region at around 50–60°N and strongly drop going northward (see Supplement Fig. S7.1). Though the diagnostic model would be able follow the spatial pattern in the climatology if data would say so (green line), in the absence of data it extrapolates the field northward keeping the high seasonal cycle from the pixels south of this area.

Seasonal data coverage is also limiting in the tropical Indian, the other problematic region (see Supplement Fig. S7.4).

## Appendix C

### Calculating phosphate concentrations from $f_{\text{int}}^{\text{DIC}}$

Considering the ocean-internal processes affecting mixed-layer concentrations of biogeochemical tracers, carbon addition/removal occurs both through biological respiration/photosynthesis in the mixed layer and through mixing-in of water masses with higher/lower carbon concentration. In the second case, in turn, there is a contribution from the pre-formed carbon concentration originating from the last contact of the arriving water parcel with the atmosphere, and a contribution from carbon added to the water parcel along its way through the ocean interior by remineralization.

The biological contribution to the sources and sinks of carbon (both in the mixed layer and the ocean interior) is directly linked to sources and sinks of phosphate, and assumed to proceed in Redfield proportions $r_{C:P} = 106$. The remainder is proportional to the differences between the pre-formed carbon (or nutrient) concentrations of arriving water parcels and those of the destination mixed-layer. Considering only their variability as relevant here, these differences are assumed to be either small or close to Redfield proportions as well (this is equivalent to assuming that the internal sources/sinks of the tracer $C_m^{\text{DIC}*} = C_m^{\text{DIC}} - r_{C:P}C_m^{\text{PO}_4}$ (Gloor et al., 2003) have no variability). Then the variations of ocean-internal phosphate sources and sinks ($f_{\text{int}}^{\text{PO}_4}$) are proportional to the variations of $f_{\text{int}}^{\text{DIC}}$. We thus calculated the potential seasonal changes in the mixed-layer $\text{PO}_4$ concentration using a $\text{PO}_4$ budget equation as Eq. (A18) (including history and freshwater fluxes), except that (1) the ocean-internal flux was reduced by the Redfield ratio ($f_{\text{int}}^{\text{PO}_4} = f_{\text{int}}^{\text{DIC}}/r_{C:P}$), (2) the long-term mean was subtracted from $f_{\text{int}}^{\text{PO}_4}$, and (3) there was no sea-air flux of phosphate ($f_{\text{ma}}^{\text{PO}_4} = 0$). As the initial condition for the budget equation is unknown and phosphate does not relax towards an atmospheric value, calculated $C_m^{\text{PO}_4}$ contains an ar-
bitrary offset, such that we only compare the deviation from
the temporal mean.

**Supplementary material related to this article is**
**available online at:** http://www.ocean-sci.net/9/193/2013/os-9-193-2013-supplement.pdf.

*Acknowledgements.* We would like to thank all contributors to the SOCAT database, as well as all providers of atmospheric $CO_2$ data, which are the basis of this work. We gratefully acknowledge helpful discussions with Nicolas Cassar, Andrew Dickson, Niki Gruber, Roberta Hamme, Steve Jones, Bob Key, Armin Köhl, Wolfgang Koeve, Iris Kriest, Corinne Le Quéré, Miguel Mahecha, Bill Munger, Andreas Oschlies, Taro Takahashi, Yasunori Tohjima, and Lisan Yu. Sara Mikaloff Fletcher and Rik Wanninkhof made very helpful comments as referees. This study was also supported through EU FP7 project CARBOCHANGE (264879) and is a contribution to the international IMBER/SOLAS projects. The service charges for this open access publication have been covered by the Max Planck Society.

Edited by: M. Hoppema

## References

Archer, D., Kheshgi, H., and E. M.-R.: Multiple timescales for neutralization of fossil fuel $CO_2$, Geophys. Res. Lett., 24, 405–408, 1997.

Arthun, M., Bellerby, R. G. J., Omar, A. M., and Schrum, C.: Spatiotemporal variability of air-sea $CO_2$ fluxes in the Barents Sea, as determined from empirical relationships and modeled hydrography, J. Mar. Sys., 98–99, 40–50, doi:10.1016/j.jmarsys.2012.03.005, 2012.

Atlas, R., Hoffman, R. N., Ardizzone, J., Leidner, S. M., Jusem, J. C., Smith, D. K., and Gombos, D.: A cross-calibrated, multi-platform ocean surface wind velocity product for meteorological and oceanographic applications, B. Am. Meteorol. Soc., 92, 157–174, 2011.

Baker, D. F., Law, R. M., Gurney, K. R., Rayner, P., Peylin, P., Denning, A. S., Bousquet, P., Bruhwiler, L., Chen, Y.-H., Ciais, P., Fung, I. Y., Heimann, M., John, J., Maki, T., Maksyutov, S., Masarie, K., Prather, M., Pak, B., Taguchi, S., and Zhu, Z.: TransCom 3 inversion intercomparison: Impact of transport model errors on the interannual variability of regional $CO_2$ fluxes, 1988–2003, Global Biogeochem. Cy., 20, GB1002, doi:10.1029/2004GB002439, 2006.

Bates, N. R., Knap, A. H., and Michaels, A. F.: Contribution of hurricanes to local and global estimates of air-sea exchange of $CO_2$, Nature, 395, 58–61, 1997.

Bousquet, P., Peylin, P., Ciais, P., Le Quéré, C., Friedlingstein, P., and Tans, P. P.: Regional changes in carbon dioxide fluxes of land and oceans since 1980, Science, 290, 1342–1346, 2000.

Buitenhuis, E., Rivkin, R., Sailley, S., and Le Quéré, C.: Biogeochemical fluxes through microzooplankton, Global Biogeochem. Cy., 24, GB4015, doi:10.1029/2009GB003601, 2010.

Chen, L., Xu, S., Gao, Z., Chen, H., Zhang, Y., Zhan, J., and Li, W.: Estimation of monthly air-sea $CO_2$ flux in the southern Atlantic and Indian Ocean using in-situ and remotely sensed data, Remote Sens. Environ., 115, 1935–1941, 2011.

Ciais, P., Tans, P. P., Trolier, M., White, J. W. C., and Francey, R. J.: A Large Northern Hemisphere Terrestrial $CO_2$ Sink Indicated by the $^{13}$C/$^{12}$C Ratio of Atmospheric $CO_2$, Science, 269, 1098–1102, doi:10.1126/science.269.5227.1098, 1995.

Conkright, M., Locarnini, R. A., Garcia, H., O’Brien, T., Boyer, T., Stephens, C., and Antonov, J.: World Ocean Atlas 2001: Objective Analyses, Data Statistics, and Figures, CD-ROM Documentation, National Oceanographic Data Center, Silver Spring, MD, 17 pp., 2002.

Conway, T., Tans, P., Waterman, L., Thoning, K., Kitzis, D., Masarie, K., and Zhang, N.: Evidence for interannual variability of the carbon cycle from the national oceanic and atmospheric administration climate monitoring and diagnostics laboratory global air sampling network, J. Geophys. Res., 99, 22831–22855, 1994.

Corbière, A., Metzl, N., Reverdin, G., Brunet, C., and Takahashi, T.: Interannual and decadal variability of the oceanic carbon sink in the North Atlantic subpolar gyre, Tellus B, 59, 168–178, 2007.

de Boyer Montégut, C., Madec, G., Fischer, A. S., Lazar, A., and Iudicone, D.: Mixed layer depth over the global ocean: an examination of profile data and a profile-based climatology, JGR, 109, C12003, doi:10.1029/2004JC002378, 2004.

Egleston, E. S., Sabine, C. L., and Morel, F. M. M.: Revelle revisited: Buffer factors that quantify the response of ocean chemistry to changes in DIC and alkalinity, Global Biogeochem. Cy., 24, GB1002, doi:10.1029/2008GB003407, 2010.

Garcia, H. E., Locarnini, R. A., Boyer, T. P., and Antonov, J. I.: World Ocean Atlas 2005, Volume 4: Nutrients (phosphate, nitrate, silicate), in: NOAA Atlas NESDIS 64, edited by Levitus, S., US Government Printing Office, Washington, DC, 2006.

Gloor, M., Gruber, N., Sarmiento, J., Sabine, C. L., Feely, R. A., and Rödenbeck, C.: A first estimate of present and preindustrial air-sea $CO_2$ flux patterns based on ocean interior carbon measurements and models, Geophys. Res. Lett., 30, 1010, doi:10.1029/2002GL015594, 2003.

Heimann, M. and Körner, S.: The global atmospheric tracer model TM3, Tech. Rep. 5, Max Planck Institute for Biogeochemistry, Jena, Germany, 2003.

Jacobson, A. R., Mikaloff Fletcher, S. E., Gruber, N., Sarmiento, J. L., and Gloor, M.: A joint atmosphere-ocean inversion for surface fluxes of carbon dioxide: 1. Methods and global-scale fluxes, Global Biogeochem. Cy., 21, GB1019, doi:10.1029/2005GB002556, 2007.

Jähne, B., Heinz, G., and Dietrich, W.: Measurement of the diffusion coefficients of sparingly soluble gases in water with a modified Barrer method, J. Geophys. Res., 92, 10767–10776, 1987.

Jones, S. D., Le Quéré, C., and Rödenbeck, C.: Autocorrelation characteristics of surface ocean $p^{CO_2}$ and air-sea $CO_2$ fluxes, Glob. Biogeochem. Cy., 26, GB2042, doi:10.1029/2010GB004017, 2012a.

Jones, S. D., Le Quéré, C., Rödenbeck, C., and Manning, A.: A Global Interpolation of Surface Ocean $p^{CO_2}$, J. Adv. Model. Earth Sys., (in review), 2012b.

Joes, F., Plattner, G.-K., Stocker, T. F., Marchal, O., and Schmittner, A.: Global warming and marine carbon cycle feedbacks an future
atmospheric CO₂, Science, 284, 464–467, 1999.

Kalnay, E., Kanamitsu, M., Kistler, R., Collins, W., Deaven, D., Gandin, L., Iredell, M., Saha, S., White, G., Woollen, J., Zhu, Y., Chelliah, M., Ebisuzaki, W., Higgins, W., Janowiak, J., Mo, K. C., Ropelewski, C., Wang, J., Leetmaa, A., Reynolds, R., Jenne, R., and Joseph, D.: The NCEP/NCAR 40-year reanalysis project, Bull. Am. Met. Soc., 77, 437–471, 1996.

Keeling, R. F. and Shertz, S. R.: Seasonal and interannual variations in atmospheric oxygen and implications for the global carbon cycle, Nature, 358, 723–727, 1992.

Key, R. M., Kozyr, A., Sabine, C. L., Lee, K., Wanninkhof, R., Bullister, J. L., Feely, R. A., Millero, F. J., Mordy, C., and Peng, T.-H.: A global ocean carbon climatology: Results from Global Data Analysis Project (GLODAP), Global Biogeochem. Cy., 18, GB4031, doi:10.1029/2004GB002247, 2004.

Krakauer, N. Y., Randerson, J. T., Primeau, F. W., Gruber, N., and Menemenlis, D.: Carbon isotope evidence for the latitudinal distribution and wind speed dependence of the air-sea gas transfer velocity, Tellus B, 58, 390–417, 2006.

Law, R. M., Peters, W., Rödenbeck, C., Aulagnier, C., Baker, I., Bergmann, D. J., Bousquet, P., Brandt, J., Bruhwiler, L., Cameron-Smith, P. J., Christensen, J. H., Delage, F., Denning, A. S., Fan, S., Geels, C., Houweling, S., Imasu, R., Karstens, U., Kawa, S. R., Kleist, J., Krol, M. C., Lin, S.-J., Lokupitiya, R., Maki, T., Maksyutov, S., Niwa, Y., Onishi, R., Parazoo, N., Patra, P. K., Pieterse, G., Rivier, L., Satoh, M., Serrar, S., Taguchi, S., Takigawa, M., Vautard, R., Vermeulen, A. T., and Zhu, Z.: TransCom model simulations of hourly atmospheric CO₂: Experimental overview and diurnal cycle results for 2002, Global Biogeochem. Cy., 22, GB3009, doi:10.1029/2007GB003050, 2008.

Le Quéré, C., Rödenbeck, C., Buitenhuis, E. T., Conway, T. J., Langenfelds, R., Gomez, A., Labuschagne, C., Ramonet, M., Nakazawa, T., Metzl, N., Gillett, N., and Heimann, M.: Saturation of the Southern Ocean CO₂ sink due to recent climate change, Science, 316, 1735–1738, 2007.

Le Quéré, C., Takahashi, T., Buitenhuis, E. T., Rödenbeck, C., and Sutherland, S. C.: Impact of climate change and variability on the global oceanic sink of CO₂, Global Biogeochem. Cy., 24, GB4007, doi:10.1029/2009GB003599, 2010.

Lee, K., Tong, L. T., Millero, F. J., Sabine, C. L., Dickson, A. G., Goyet, C., Park, G.-H., Wanninkhof, R., Feely, R. A., and Key, R. M.: Global relationships of total alkalinity with salinity and temperature in surface waters of the world's oceans, Geophys. Res. Lett., 33, L19605, doi:10.1029/2006GL027207, 2006.

Lefèvre, N., Watson, A. J., and Watson, A. R.: A comparison of multiple regression and neural network techniques for mapping in situ pCO₂ data, Tellus B, 75, 375–384, 2005.

Lourantou, A. and Metzl, N.: Decadal evolution of carbon sink within a strong bloom area in the subantarctic zone, Geophys. Res. Lett., 38, L23608, doi:10.1029/2011GL049614, 2011.

Matear, R. J. and Hirst, A. C.: Climate change feedback on the future oceanic CO₂ uptake, Tellus B, 51, 722–733, 1999.

Metzl, N., Corbière, A., Reverdin, G., Lenton, A., Takahashi, T., Olsen, A., Johannessen, T., Pierrot, D., Wanninkhof, R., Ólafsdóttir, S. R., Olafsson, J., and Ramonet, M.: Recent acceleration of the sea surface fCO₂ growth rate in the North Atlantic subpolar gyre (1993–2008) revealed by winter observations, Global Biogeochem. Cy., 24, GB4004, doi:10.1029/2009GB003658, 2010.

Mikaloff Fletcher, S. E., Gruber, N., Jacobson, A. R., Doney, S. C., Dutkiewicz, S., Gerber, M., Follows, M., Joos, F., Lindsay, K., Menemenlis, D., Mouchet, A., Müller, S. A., and Sarmiento, J. L.: Inverse estimates of anthropogenic CO₂ uptake, transport, and storage by the ocean, Global Biogeochem. Cy., 20, GB2002, doi:10.1029/2005GB002530, 2006.

Müller, S. A., Joos, F., Plattner, G.-K., Edwards, N. R., and Stocker, T. F.: Modelled natural and excess radiocarbon - sensitivities to the gas exchange formulation and ocean transport strength, Global Biogeochem. Cy., 22, GB3011, doi:10.1029/2007GB003065, 2008.

Naegler, T.: Reconciliation of excess ¹⁴C-constrained global CO₂ piston velocity estimates, Tellus B, 61, 372–384, 2009.

Naegler, T. and Levin, I.: Closing the global radiocarbon budget 1945–2005, J. Geophys. Res., 111, D12311, doi:10.1029/2005JD006758, 2006.

Newsam, G. N. and Enting, I. G.: Inverse problems in atmospheric constituent studies: I. Determination of surface sources under a diffusive transport approximation, Inverse Problems, 4, 1037–1054, 1988.

Olivier, J. G. J. and Berdowski, J. J. M.: Global emissions sources and sinks, in: The Climate System, Berdowski, J., Guicherit, R., and Heij, B. J., A. A. Balkema Publishers/Swets and Zeitlinger Publishers, Lisse, The Netherlands, ISBN 9058092550, 33–78, 2001.

Park, G.-H., Wanninkhof, R., Doney, S. C., Takahashi, T., Lee, K., Feely, R. A., Sabine, C. L., Triñanes, J., and Lima, I. D.: Variability of global net sea-air CO₂ fluxes over the last three decades using empirical relationships, Tellus B, 62, 352–368, 2010.

Pfeil, B., Olsen, A., Bakker, D. C. E., Hankin, S., Koyuk, H., Kozyr, A., Malczyk, J., Manke, A., Metzl, N., Sabine, C. L., Akl, J., Alin, S. R., Bellerby, R. G. J., Borges, A., Boutin, J., Brown, P. J., Cai, W.-J., Chavez, F. P., Chen, A., Cosca, C., Fassbender, A. J., Feely, R. A., González-Dávila, M., Goyet, C., Hardman-Mountford, N., Heinze, C., Hood, M., Hoppema, M., Hunt, C. W., Hydes, D., Ishii, M., Johannessen, T., Jones, S. D., Key, R. M., Körtzinger, A., Landschützer, P., Lauvset, S. K., Lefèvre, N., Lenton, A., Lourantou, A., Merlivat, L., Midorikawa, T., Mintrop, L., Miyazaki, C., Murata, A., Nakadate, A., Nakano, Y., Nakaoka, S., Nojiri, Y., Omar, A. M., Padin, X. A., Park, G.-H., Paterson, K., Perez, F. F., Pierrot, D., Poisson, A., Ríos, A. F., Santana-Casiano, J. M., Salisbury, J., Sarma, V. V. S. S., Schlitzer, R., Schneider, B., Schuster, U., Sieger, R., Skjelvan, I., Steinhoff, T., Suzuki, T., Takahashi, T., Tedesco, K., Telszewski, M., Thomas, H., Tilbrook, B., Tjiputra, J., Vande-mark, D., Veness, T., Wanninkhof, R., Watson, A. J., Weiss, R., Wong, C. S., and Yoshikawa-Inoue, H.: A uniform, quality controlled Surface Ocean CO₂ Atlas (SOCAT), Earth Syst. Sci. Data Discuss., 5, 735–780, doi:10.5194/essdd-5-735-2012, 2012.

Rayner, P., Enting, I., Francey, R., and Langenfelds, R.: Reconstructing the recent carbon cycle from atmospheric CO₂, δ¹³CO₂ and O₂/N₂ observations, Tellus B, 51, 213–232, 1999.

Rödenbeck, C.: Estimating CO₂ sources and sinks from atmospheric mixing ratio measurements using a global inversion of atmospheric transport, Tech. Rep. 6, Max Planck Institute for Biogeochemistry, Jena, Germany, 2005.

Rödenbeck, C., Houweling, S., Gloor, M., and Heimann, M.: CO₂ flux history 1982–2001 inferred from atmospheric data using a
global inversion of atmospheric transport, Atmos. Chem. Phys.,
3, 1919–1964, doi:10.5194/acp-3-1919-2003, 2003.

Sabine, C., Maenner, S. and Sutton, A.: High-resolution ocean and atmosphere pCO2 time-series measurements from mooring TAO140W. http://cdiac.esd.ornl.gov/ftp/oceans/Moorings/TAO140W/. Carbon Dioxide Information Analysis Center, Oak Ridge National Laboratory, US Department of Energy, Oak Ridge, Tennessee, 2010.

Sarmiento, J. L. and Gruber, N.: Ocean Biogeochemical Dynamics, Princeton Univ. Press, Princeton, USA, 2006.

Sarmiento, J. L. and Le Quéré, C.: Oceanic carbon dioxide uptake in a model of century-scale global warming, Science, 274, 1346–1350, 1996.

Sarmiento, J. L., Gloor, E., Gruber, N., Beaulieu, C., Jacobson, A. R., Fletcher, S. E. M., Pacala, S., and Rodgers, K.: Trends and Regional Distributions of Land and Ocean Carbon Sinks, Biogeosciences, 7, 2351–2367, doi:10.5194/bg-7-2351-2010, 2010.

Sigman, D. M. and Boyle, E. A.: Glacial/interglacial variations in atmospheric carbon dioxide, Nature, 407, 859–869, 2000.

Stephens, B. B., Gurney, K. R., Tans, P. P., Sweeney, C., Peters, W., Bruhwiler, L., Ciais, P., Ramonet, M., Bousquet, P., Nakazawa, T., Aoki, S., Machida, T., Inoue, G., Vinnichenko, N., Lloyd, J., Jordan, A., Heimann, M., Shibistova, O., Langenfelds, R. L., Steele, L. P., Francey, R. J., and Denning, A. S.: Weak northern and strong tropical land carbon uptake from vertical profiles of atmospheric CO₂, Science, 316, 1732–1735, 2007.

Sweeney, C., Gloor, E., Jacobson, A. R., Key, R. M., McKinley, G., and co authors: Constraining global air-sea gas exchange for CO₂ with recent Bomb ¹⁴C measurements, Global Biogeoch. Cycl., 21, GB2015, doi:10.1029/2006GB002784, 2007.

Takahashi, T., Sutherland, S. C., Wanninkhof, R., Sweeney, C., Feely, R. A., Chipman, D. W., Hales, B., Friederich, G., Chavez, F., Sabine, C., Watson, A., Bakker, D. C. E., Schuster, U., Metzl, N., Yoshikawa-Inoue, H., Ishii, M., Midorikawa, T., Nojiri, Y., Kortzinger, A., Steinhoff, T., Hoppema, M., Olafsson, J., Arnarson, T. S., Tillbrook, B., Johannessen, T. and Olsen, A., Bellerby, R., Wong, C. S., Delille, B., Bates, N. R., and de Baar, H. J. W.: Climatological mean and decadal change in surface ocean pCO₂ and net sea-air CO₂ flux over the global oceans, Deep-Sea Res. II, 56, 554–577, 2009.

Takahashi, T., Sutherland, S. C., and Kozyr, A.: Global Ocean Surface Water Partial Pressure of CO₂ Database: Measurements Performed during 1968–2010 (Version 2010), Tech. Rep. NDP-088, ORNL/CDIAC-152, 2010.

Telszewski, M., Chazottes, A., Schuster, U., Watson, A., Moulin, C., Bakker, D., González-Dávila, M., Johannessen, T., Körtzinger, A., Lüger, H., Olsen, A., Omar, A., Padin, X., Ríos, A., Steinhoff, T., Santana-Casiano, M., Wallace, D. W. R., and Wanninkhof, R. H.: Estimating the monthly pCO₂ distribution in the North Atlantic using a self-organizing neural network, Biogeosciences, 6, 1405–1421, doi:10.5194/bg-6-1405-2009, 2009.

Valsala, K. V. and Maksyutov, S.: Simulation and assimilation of global ocean pCO₂ and air-sea CO₂ fluxes using ship observations of surface ocean pCO₂ in a simplified Biogeochemical offline model, Tellus B, 62, 821–840, doi:10.1111/j.1600-0889.2010.00495.x, 2010.

Wanninkhof, R.: Relationship between wind speed and gas exchange over the ocean, J. Geophys. Res.-Oceans, 97, 7373–7382, 1992.

Wanninkhof, R., Park, G.-H., Takahashi, T., Sweeney, C., Feely, R., Nojiri, Y., Gruber, N., Doney, S. C., McKinley, G. A., Lenton, A., Le Quéré, C., Heinze, C., Schwinger, J., Graven, H., and Khatiwala, S.: Global ocean carbon uptake: magnitude, variability and trends, Biogeosciences Discuss., 9, 10961–11012, doi:10.5194/bgd-9-10961-2012, 2012.

Weiss, R.: Carbon dioxide in water and seawater: the solubility of a non-ideal gas., Mar. Chem., 2, 203–205, 1974.

Weiss, R. F.: The solubility of nitrogen, oxygen and argon in water and seawater, Deep-Sea Res., 17, 721–735, 1970.

Weiss, R. F. and Price, B. A.: Nitrous oxide solubility in water and seawater, Mar. Chem., 8, 347–359, 1980.

While, J., Totterdell, I., and Martin, M.: Assimilation of pCO₂ Data Into a Global Coupled Physical-Biogeochemical Ocean Model, J. Geophys. Res., 117, C03037, doi:10.1029/2010JC006815, 2012.

Yu, L. and Weller, R. A.: Objectively Analyzed air-sea heat Fluxes (OAFlux) for the global oceans, Bull. Am. Met. Soc., 88, 527–539, 2007.

Zeebe and Wolf-Gladrow: CO₂ in Seawater: Equilibrium, Kinetics, Isotopes, Elsevier, Amsterdam, 2001.