Global surface-ocean p CO 2 and sea – air CO 2 flux variability from an observation-driven ocean mixed-layer scheme

A temporally and spatially resolved estimate of the global surface-ocean CO 2 partial pressure field and the sea–air CO2 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 CO2 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.


Introduction
The oceans are considered the dominant player in the global carbon cycle on long timescales, e.g. in the glacialinterglacial 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. 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: 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 2 exchange between the Earth surface and the atmosphere can be quantified based on atmospheric CO 2 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 2 exchange is much less variable and smoothed out in time because the carbonate chemistry of seawater slows the equilibration rate of dissolved CO 2 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 2 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 seaair CO 2 fluxes over large spatial regions, not its temporal or high-resolution spatial variability. The long-term global seaair CO 2 flux has also been estimated from observed trends in atmospheric oxygen (Keeling and Shertz, 1992) as well as from 13 C isotopic ratios in atmospheric CO 2 (Ciais et al., 1995).
In view of the above-mentioned problem of the atmospheric CO 2 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 databased sea-air CO 2 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 2 Atlas (SOCAT) database (version 1.5) of CO 2 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 2 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 4 ) climatology. Estimated interannual variations of the sea-air CO 2 flux are presented and discussed in a companion paper (Rödenbeck et al., 2013).

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 2 exchange (f CO 2 ma ) 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 f DIC int (x,y,t) 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.
Atm. Inv. Fig. 2. Illustration of the inverse procedure (three modi). Boxes denote the process parameterizations from Figure 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 CO 2 ma are adjusted to match atmospheric mixing ratios (terrestrial Net Ecosystem Exchange is adjusted as well but omitted from this schematic for clarity). Middle: In case ATM, ocean-internal sources/sinks f DIC int are adjusted instead of sea-air fluxes to match atmospheric CO2 mixing ratios (again Net Ecosystem Exchange is further adjusted). Right: In case SFC, ocean-internal sources/sinks are adjusted to match the surface-ocean CO2 partial pressure observations. Sea-air fluxes are then calculated from the estimated partial pressure field. 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. 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 seasurface CO 2 partial pressure. Through cost function contributions gauging the match of modelled and measured p CO 2 , 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).

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 2 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-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.

Atm. Inv.
For ATM For SFC Fig. 2. Illustration of the inverse procedure (three modi). Boxes denote the process parameterizations from Figure 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 CO 2 ma are adjusted to match atmospheric mixing ratios (terrestrial Net Ecosystem Exchange is adjusted as well but omitted from this schematic for clarity). Middle: In case ATM, ocean-internal sources/sinks f DIC int are adjusted instead of sea-air fluxes to match atmospheric CO2 mixing ratios (again Net Ecosystem Exchange is further adjusted). Right: In case SFC, ocean-internal sources/sinks are adjusted to match the surface-ocean CO2 partial pressure observations. Sea-air fluxes are then calculated from the estimated partial pressure field.

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 CO 2 ma 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 DIC int are adjusted instead of sea-air fluxes to match atmospheric CO 2 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 2 partial pressure observations. Sea-air fluxes are then calculated from the estimated partial pressure field. 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 2 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 2 ) only account for a small part of the carbon relevant in the ocean-internal budget (dissolved inorganic carbon, DIC). The link between CO 2 abundance (expressed in terms of partial pressure p CO 2 m ) and DIC abundance (in terms of its concentration C DIC m ) is determined by chemical equilibria, which we assume to be attained instantaneously. The non-linear dependence of p CO 2 m on C DIC m has been linearized in standard ways . The chemistry parameterization also contains a temperature-dependent factor, and contributions from seasonal variations in alkalinity and salinity (Appendix A1.2).
www.ocean-sci.net/9/193/2013/ Ocean Sci., 9, 193-216, 2013 Mixed-layer carbon budget. Changes in the spatiotemporal 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 firstorder differential equation. Its solution gives DIC concentration C DIC m as a function of the ocean-internal sinks f DIC int (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.

Main adjustable degrees of freedoms
The ocean-internal carbon sources and sinks (f DIC int ) 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.

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.

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).

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 oceaninterior 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.

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 DIC int, pri = 0 (Appendix A2.2, see bottom panel). The choice f DIC int, pri = 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 DIC int . The variations in the prior values of C DIC m , p CO 2 , and the sea-air fluxes f CO 2 ma , 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 DIC int 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 1): Sea-air CO2 fluxes (f CO 2 ma ), surface-ocean CO2 partial pressure (p CO 2 m ), mixed-layer DIC concentration (C DIC m ), as well as ocean-internal carbon addition to (removal from) the mixed layer (f DIC int ). Estimates have been obtained by fitting the diagnostic scheme to SOCAT surface-ocean CO2 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 CO 2 m panel. The Bayesian prior (f DIC int = 0 and corresponding C DIC m , p CO 2 m , and f CO 2 ma , Appendix A2.2, Sect. 3.2) is also given (thin dashed gray).  , as well as ocean-internal carbon addition to (removal from) the mixed layer (f DIC int ). 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   Figure 1): Sea-air CO 2 fluxes (f CO 2 ma ), surface-ocean CO 2 partial pressure (p CO 2 m ) ocean-internal carbon addition to (removal from) the mixed layer (f DIC int ). Estimates hav SOCAT surface-ocean CO 2 partial pressure data (SFC, blue). The SOCAT data points th black dots in the p CO 2 m panel. The Bayesian prior (f DIC int = 0 and corresponding C DIC m , p CO 2 m (thin dashed gray).

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 CO 2 m at each pixel. Constrained Extrapolated pixel pixel Where data values exist, pixels are directly constrained (left column, which is identical to Figure 2). Due to the a-priori spatial correlations (Appendix A2.2), adjustments to f DIC int at the constrained pixel also change f DIC int at neighbouring unconstrained pixels. p CO 2 is then calculated passively from the extrapolated f DIC int by the process parameterizations, using the local values of the driving fields (right column). In data-dense areas, local adjustments to f DIC int (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, Figure 13). 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 DIC int ). In data-poor areas, there may exist pixels that do not have any constrained pixel within the correlation radius; then f DIC int 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). 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 DIC int at the constrained pixel also change f DIC int at neighbouring unconstrained pixels. p CO 2 is then calculated passively from the extrapolated f DIC int by the process parameterizations, using the local values of the driving fields (right column). In data-dense areas, local adjustments to f DIC int (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 DIC int ). In data-poor areas, there may exist pixels that do not have any constrained pixel within the correlation radius; then f DIC int 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 apriori correlations defined in Appendix A2.2. The mechanism is illustrated in Fig. 5. In particular, as seasonal variations in f DIC int are allowed to be adjusted much more freely than non-seasonal variations (equivalent to a-priori correlations between adjustments to f DIC int 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 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 interannual variations requires more even data distribution in time than constraining seasonal variability; this density of data is only available is some areas of the ocean (Rödenbeck et al., 2013).
The spatio-temporal a-priori correlations also act to smooth f DIC int 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 DIC int the corresponding p 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 CO 2 fields from a biogeochemical process model simulation (as the model has less day-to-day variability than the real p CO 2 field, however, this test may underestimate errors from aliasing such high-frequency variations into variations on the resolved timescales).

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 DIC int 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 CO 2 . In most regions, the sensitivity of the sea-air flux is slightly larger than that of p CO 2 , almost entirely due to the uncertainty of gas exchange (variation in piston velocity). Larger spread in p CO 2 than in the sea-air flux 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).

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  from scale mismatch between the model pixel (≈ 4 • lat. × 5 • 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 • -90 • 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.   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 shortterm 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 CO 2 m as reproduced by the solubility and chemistry parameterizations, tropical p CO 2 m 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 DIC int , which however cannot be constrained from the SO-CAT 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 Ocean Sci., 9, 193-216, 2013 www.ocean-sci.net/9/193/2013/ chemistry. In the tropics, more realistic short-term variations could potentially be obtained by parameterization of fast processes (upwelling events) in f DIC int .

Comparison to the Takahashi climatology
The monthly mean seasonal cycle calculated from our spatiotemporal 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 DIC int 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 equa-tion 1 . In data-dense areas, results should not depend much on the extrapolation method.
In most regions, the seasonal cycle estimated from SO-CAT 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 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.
www.ocean-sci.net/9/193/2013/ Ocean Sci., 9, 193-216, 2013 Table 2. Mathematical symbols (also see Table 3). Some quantities are indexed by species or other distinctions explained in the text. 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).

Quantity Unit a Meaning
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).  2010)). Comparison between the measurements (black dots) and the p CO 2 field estimated by fitting to the SOCAT data (blue). The estimated field has been picked at the times where measurements exist; one year as been selected for each site. Bottom: Sea surface temperature measured simultaneously with p CO 2 (black dots) compared to our SST driving field at this pixel (light grey).  . Comparison between the measurements (black dots) and the p CO 2 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 p CO 2 (black dots) compared to our SST driving field at this pixel (light grey).

Ocean
Beyond seasonality, agreement is also found in the im-  Table 2, "Net flux" and "Coastal area") using CCMP winds (Atlas et al., 2011).
In summary, we take the similarity between the seasonalities from the two methods based on similar data as an indication that the estimated p CO 2 field in most regions is informed by the data and not primarily method-dependent. Figure 11 compares the SOCAT-based fit from Fig. 6 (run SFC) with the fit of the diagnostic scheme to atmospheric CO 2 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 uncertain model details, not shown) than those estimated from the p CO 2 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 influence, while the large disagreement in the temperate North Pacific is due to the short but strong sink period during summer, more typical for terrestrial boreal ecosystems.

Comparison of p CO 2 -based and atmospheric results
The results establish that oceanic and atmospheric data are consistent in regions away from land influences, but that estimates based on p CO 2 data provide much stronger constraints on internal ocean processes, in particular in ocean areas closer to land. Combining the oceanic and atmospheric constraints is therefore both possible and beneficial (compare Jacobson et al., 2007). It can be implemented by adding together the cost function contributions of p CO 2 data (from SFC) and of atmospheric CO 2 (from ATM). The estimated fields in this combined case (not shown) are almost identical to the SFC run, as expected from the weakness of the atmospheric constraint on oceanic fluxes. However, as the atmospheric 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 using the sea-air fluxes of the SFC run as priors in a subsequent "classical" atmospheric inversion (similar as in the joint inversion by Jacobson et al., 2007). Details will be discussed further in the context of interannual variations in the companion paper (Rödenbeck et al., 2013). 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 PO4 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 Figure 7; note however that calculated PO4 is also sensitive to several further model elements not included in this range. Fig. 11. Monthly mean seasonal cycle of sea-air CO 2 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 2 mixing ratio data (run ATM, magenta), as well as to the prior (no data constraint, thin dashed grey)

Plausibility of the ocean-internal carbon sources and sinks
As the diagnostic scheme (in run SFC) is constructed to match the CO 2 partial pressure data, any errors in the parameterizations linking p CO 2 and f DIC int (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 p CO 2 (and sea-air flux) fields will be only little affected by imperfections in these parameterizations. Even at pixels not directly constrained by a p CO 2 data value but only indirectly via extrapolation of the f DIC int 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 p CO 2 estimate shown in Sect. 3.3. Therefore, as 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 PO4 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 Figure 7; note however that calculated PO4 is also sensitive to several further model elements not included in this range.

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 4 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 4 is also sensitive to several further model elements not included in this range.
long as we are only interested in the p CO 2 (and sea-air flux) fields, f DIC int could be just considered a mathematical device in the mapping of p CO 2 .
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 DIC int ), correctness of the parameterizations is needed. On the seasonal timescale considered here, the plausibility of the estimated f DIC int field can be tested in the light of mixedlayer phosphate concentrations. As explained in Appendix C, a modelled mixed-layer phosphate field can approximately be calculated from our f DIC int estimate, and can be compared to an observation-based monthly climatology (World Ocean Atlas, Garcia et al., 2006). The inferred PO 4 concentration Ocean Sci., 9, 193-216, 2013 www.ocean-sci.net/9/193/2013/ seasonality is broadly similar in phase and amplitude to the observations (Fig. 12), including characteristic region-toregion 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).

The diagnostic mixed-layer scheme in the context of other p CO 2 m mapping methods
The suitability of individual p CO 2 m 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 datadominated way, with the least complex model possible.
Statistical methods like the reference p CO 2 m 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 CO 2 m 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 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 CO 2 m 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 welldefined 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 timeindependent 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.

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. 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 biogeochemistry (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.

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, (atmospheric sign convention: positive = source to atmosphere), where p CO 2 a is the CO 2 partial pressure in the atmosphere, and p CO 2 m the CO 2 partial pressure in the homogeneous mixed layer 3 . The temperature-dependent solubility L CO 2 relates the partial pressure difference to a dissolved CO 2 concentration difference. The seawater density 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 CO 2 . We use the formulation according to Wanninkhof (1992): 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 is chosen such that the global mean CO 2 piston velocity k CO 2 Glob 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 CO 2 a = X CO 2 p dry (A3) inherits variability not only from atmospheric dry-air pressure p dry (Supplement Sect. S1.1) but also from the atmospheric CO 2 molar mixing ratio X 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 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 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 CO 2 field is small (see Supplement Sect. S2). 4 Flux pulse for atmospheric initial condition (run ATM only) Zero 3 lat. bands -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 2 partial pressure in the mixed layer is a function of mixed-layer DIC concentration C DIC m , alkalinity A, temperature T , and salinity S: In the diagnostic scheme, the relevant dependence is that on C DIC m , while the influences of A, T , and S are calculated apriori 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 CO 2 m at a reference temperature T Ref : For iso-chemical sea water (i.e. constant C DIC m , A, and S), 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 CO 2 m are non-linear, but monotonic and can be linearized to good approximation . Linearization of the dependence on C DIC m is important here in order to be able to use the fast minimization algorithm of Rödenbeck (2005) (Sect. A2.1). We consider deviations ( ) of the driver variables from temporally constant reference values (superscript "Ref", see Sect. A1.4 below): Linearization of Eq. (A5) around these references gives The DIC sensitivity has been calculated from the response factors γ DIC by Egleston et al. (2010): The alkalinity sensitivity has been calculated using the approximation (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 Sarmiento and Gruber, 2006). Any seasonal and interannual variations in the proportionalities ∂p CO 2 m /∂C DIC m , etc., are neglected. The variations in alkalinity, A, and salinity, S, are fixed according to input data sets (Table 1). Their effect on p CO 2 m 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 CO 2 m and C DIC m 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, Using this definition, Eq. (A10) becomes Now also defining an "apparent piston velocity" of DIC as Eqs. (A1) and (A15) can be combined and simplified into 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: Fluxes (carbon exchange per unit horizontal area and unit time) comprise loss through sea-air exchange f CO 2 ma , seasonal re-entrainment f DIC hist and freshwater effects f DIC frw explained below, and fluxes f DIC int 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 of surface water (assumed constant, Table 3) is needed as C DIC m 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 DIC int but are parameterized explicitly 5 in the budget (Eq. A18): -The freshwater effect f DIC frw describes the dilution or enhancement of the mixed-layer DIC concentration by precipitation, evaporation, river freshwater input, or seaice formation/melting. Analogously to the often used salinity normalization , 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: This will tend to overestimate the effect to the extent that freshwater input does carry carbon.
-The "history" flux f DIC hist 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 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 DIC int during the fit to the data. However, splitting them off allows easier interpretation of f DIC int (Sect. 3.7), and also prepares for the envisaged link to oxygen.
Ocean Sci., 9, 193-216, 2013 www.ocean-sci.net/9/193/2013/ meanwhile, this re-entrainment is equivalent to a carbon flux which is proportional to the difference between the current concentration C DIC m (t) and that at time t prev when the mixed layer was as deep for the last time. We call this the "history" flux. It can be written as The Heaviside step function ensures that f 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 DIC int with the mean sea-air exchange f CO 2 ma (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 DIC int . Substituting Eq. (A17) into Eq. (A18) gives where we further exploited that the temporally constant reference C DIC, Ref m disappears from the left-hand side. The reference also cancels out in the history flux (Eq. A20), which is thus a linear functional of C DIC m (t). Equation (A21) thus represents a 1 st order differential equation in C DIC m that can be solved (numerically) for given "forcing" f DIC int . At the beginning t i of the time period, we use the initial condition

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 2 inversion") described in Rödenbeck (2005) and reviewed in Sect. A2.1. This atmospheric transport inversion is extended by representing seaair fluxes by the ocean parameterizations (Sect. A2.2) and by using the p 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 obs (time series at observation sites, see Table 1). The inversion calculation seeks fluxes f that lead to the best match between c obs and corresponding modelled mixing ratios c mod (f ), in the sense that the value of the cost function 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 with the transport matrix 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 www.ocean-sci.net/9/193/2013/ Ocean Sci., 9, 193-216, 2013 estimation. They are implemented here by writing the fluxes as a function of a set of parameters x: 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 This is equivalent to assuming fluxes with a-priori mean f pri and a-priori covariance matrix Q f , pri = FF T . The specification of the flux model elements f pri and F is detailed in Rödenbeck (2005). The flux is first written as a sum of contributions from several components: (f ma = sea-air flux, f nee = terrestrial net ecosystem exchange (NEE), f fos = fossil fuel burning emissions, f ini = flux pulse at beginning of inversion period to adjust initial atmospheric mixing ratio field). Each component i is represented by a structure as Eq. (A26), i.e.
where the vectors x i denote subsets of parameters in x, and the matrices F i denote the corresponding groups of columns of F. The columns of 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 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 ma ) while suppressing flux adjustments elsewhere. The detailed flux model settings for the components f i of the surface-to-atmosphere CO 2 flux in the "Jena inversion" are f 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 nee, lt + f nee, seas + f 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 highfrequency 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 √ 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 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 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 ma is not estimated directly but implemented as a function of the ocean-internal flux f int using the parameterizations of Sect. A1. Then, instead of f ma , the ocean-internal flux f int is adjusted to match the data constraint. The detailed flux-model specifications for f int are similar to those for the sea-air flux in the atmospheric inversion: -The a-priori state has been defined by f DIC int, pri = 0 (no internal sources or sinks). The "error" of this prior model f DIC int, pri = 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 Tmol yr −1 (with respect to 3-monthly anomalies).
-Bayesian a-priori spatial and temporal correlations have been implemented that dampen variability in f DIC int on scales smaller than around 640 km (longitude), 320 km   Top: map of a-priori correlation coefficients with respect to an example pixel at (70 • E, 49 • S), illustrating how inverse adjustments to f DIC int 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) 6 , and similar to those of p CO 2 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 p CO 2 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 datasparse Southern Ocean as they lead to slightly better performance in synthetic-data tests (analogously to Appendix B, test not shown). In any case, simple distancebased 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 nonseasonal variability by a factor 16 with respect to the seasonal variability (this is accomplished by implementing the time dependence of f DIC int as a Fourier series and reducing the a-priori uncertainty of the non-seasonal www.ocean-sci.net/9/193/2013/ Ocean Sci., 9, 193-216, 2013 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 SO-CAT data in the individual years, and extrapolate the mean seasonal cycle of f DIC int 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 C DIC, ini m , needed only for technical reasons: Through adjusting C DIC, ini m , the initial conditions of the carbon budget equation can become consistent with the data (see Supplement Sect. S1.2). Adjustments to C DIC, ini m 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 CO 2 m 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 SO-CAT 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: 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 √ 2 × 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 CO 2 m -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 + 1 2 x T x).

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. 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 gray, partially off-scale) is given. The p CO 2 fields have been averaged over the regions as in Figure 7.

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  (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.
-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).

Calculating phosphate concentrations from f DIC int
Considering the ocean-internal processes affecting mixedlayer 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 preformed 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 DIC* m = C DIC m − r C:P C PO 4 m (Gloor et al., 2003) have no variability). Then the variations of ocean-internal phosphate sources and sinks (f PO 4 int ) are proportional to the variations of f DIC int . We thus calculated the potential seasonal changes in the mixed-layer PO 4 concentration using a 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 PO 4 int = f DIC int /r C:P ), (2) the long-term mean was subtracted from f PO 4 int , and (3) there was no sea-air flux of phosphate (f PO 4 ma = 0). As the initial condition for the budget equation is unknown and phosphate does not relax towards an atmospheric value, calculated C PO 4 m contains an arwww.ocean-sci.net/9/193/2013/ Ocean Sci., 9, 193-216, 2013 bitrary offset, such that we only compare the deviation from the temporal mean.