potential oxygen (APO) and land-ocean carbon sink partitioning

Abstract. A three dimensional, time-evolving field of atmospheric potential oxygen (APO ~O2/N2+CO2) was estimated using surface O2, N2 and CO2 fluxes from the WHOI ocean ecosystem model to force the MATCH atmospheric transport model. Land and fossil carbon fluxes were also run in MATCH and translated into O2 tracers using assumed O2:CO2 stoichiometries. The modeled seasonal cycles in APO agree well with the observed cycles at 13 global monitoring stations, with agreement helped by including oceanic CO2 in the APO calculation. The modeled latitudinal gradient in APO is strongly influenced by seasonal rectifier effects in atmospheric transport. An analysis of the APO-vs.-CO2 mass-balance method for partitioning land and ocean carbon sinks was performed in the controlled context of the MATCH simulation, in which the true surface carbon and oxygen fluxes were known exactly. This analysis suggests uncertainty of up to ±0.2 PgC in the inferred sinks due to variability associated with sparse atmospheric sampling. It also shows that interannual variability in oceanic O2 fluxes can cause large errors in the sink partitioning when the method is applied over short timescales. However, when decadal or longer averages are used, the variability in the oceanic O2 flux is relatively small, allowing carbon sinks to be partitioned to within a standard deviation of 0.1 Pg C/yr of the true values, provided one has an accurate estimate of long-term mean O2 outgassing.


Introduction
Atmospheric O 2 /N 2 measurements provide complementary information about atmospheric CO 2 due to the close coupling between oxygen and carbon fluxes during terrestrial photosynthesis and respiration (Keeling et al., 1993).In contrast, oceanic O 2 and CO 2 fluxes are more or less decoupled, for reasons described below.The tracer atmospheric potential oxygen (APO) exploits these differences to remove the terrestrial contribution to changes in atmospheric O 2 /N 2 , thus allowing the oceanic contribution to be largely isolated.APO is defined as, where 1.1 derives from the approximate -O 2 :CO 2 molar ratio of terrestrial photosynthesis and respiration (Severinghaus, 1995;Stephens et al., 1998;Keeling et al., 1998).The parameter X O 2 =0.2098 is the mole fraction of O 2 in dry air and converts from ppm O 2 to per meg, the units used for O 2 /N 2 (Keeling and Shertz, 1992).Oxygen and carbon fluxes are also closely coupled during fossil fuel combustion, but with a slightly larger -O 2 :CO 2 molar ratio of ∼1.4.As a result, fossil fuel combustion exerts a small influence on APO, since it yields changes in O 2 /N 2 and CO 2 that cancel out mostly but not completely in Eq. 1.
APO data have been used in at least three important carbon cycle applications.First, the seasonal cycle in APO reflects seasonal imbalances between photosynthetic new production in the surface ocean, which occurs largely in spring and summer, and the wintertime ventilation of deeper waters, which are depleted in O 2 due to remineralization.The amplitude of the seasonal cycle in APO has been used to infer the rate of seasonal new production at the ocean basin or hemispheric scale (Bender et al., 1996;Balkanski et al., 1999;Najjar and Published by Copernicus Publications on behalf of the European Geosciences Union.C. D. Nevison et al.: Air-sea O 2 flux variability and its impact on APO Keeling, 2000), although temporal and spatial overlap between new production and ventilation cause uncertainties in the inferred rates (Nevison et al., 2005;Jin et al., 2007).
Second, the latitudinal gradient in APO in principle reflects global-scale patterns of ocean biogeochemistry and ocean circulation, with O 2 uptake occurring primarily at high latitudes and outgassing of both O 2 and CO 2 taking place in the tropics.The ocean-induced patterns are superimposed on the overall north-to-south increase in APO due to northern hemisphere-dominated fossil fuel combustion.Early atmospheric transport model simulations of the latitudinal gradient in APO found a mismatch with observations, and suggested that the cause might be deficiencies in the physics of the ocean carbon models used to force the simulations (Stephens et al., 1998).Later studies suggested that the discrepancies might be attributed mainly to the atmospheric transport models (Gruber et al., 2001;Naegler et al., 2007).A further complication is that the observed APO gradient appears to be evolving with time and displays significant interannual variability (Battle et al., 2006).
A third, critically important application of APO data is to constrain the partitioning of anthropogenic CO 2 uptake between the ocean and the land biosphere, which together have absorbed more than half of the anthropogenic carbon put in the atmosphere over the last half century (Battle et al., 2000;Bender et al., 2005;Manning and Keeling, 2006; referred to hereafter as MK06).This application involves solving a system with two equations and two unknowns, as derived from the mass balances for CO 2 and O 2 /N 2 .
CO 2 / t = β(F fuel − F land − F ocean ) (2) where β=0.471 is the conversion from Pg C to ppm CO 2 , γ =1/X O 2 is the ppm to per meg conversion factor described above, Z eff is an oceanic O 2 outgassing term in per meg units, and α bio =1.1 and α f =1.4 are the -O 2 :CO 2 molar ratios of terrestrial respiration and photosynthesis and fossil fuel combustion.F fuel is the release of fossil carbon to the atmosphere, which is known from industry data (Marland et al., 2002), while F land and F ocean , the two unknowns, are CO 2 uptake by land and ocean, defined here as positive when they act as sinks for atmospheric CO 2 .Eqs. 2 and 3 can be combined according to the definition of APO in Eq. 1 to yield the APO time derivative: One can solve for F land and F ocean using either the system of equations 2 and 3 or 2 and 4 based on the global rates of change of CO 2 and O 2 /N 2 or CO 2 and APO, respectively, observed by global atmospheric monitoring networks.
The above method, originally presented as a vector diagram (Keeling et al., 1996), assumes that the O 2 :CO 2 stoichiometries of fossil fuel combustion and terrestrial respiration and photosynthesis are known, but that oceanic O 2 and CO 2 fluxes are decoupled.The vector method is often considered the best way to monitor the partitioning of sinks for anthropogenic CO 2 on a continual basis, allowing the detection of trends or changes, e.g., the increased land sink of the 1990s relative to the 1980s (Prentice et al., 2001).Alternative methods for estimating ocean and land sinks are generally more difficult to update quickly.Ocean inventory methods, for example, require a large new set of ocean cruise data (e.g., Matsumoto and Gruber, 2005).
Despite the advantages of the APO vs. CO 2 method, its inherent uncertainties may be so large that it cannot constrain the carbon sink partitioning to within better than ∼±0.7 PgC/yr (LeQuere et al., 2003).This uncertainty represents a substantial fraction of the mean total land and ocean sinks, which are estimated at 1.9 and 1.2 PgC/yr, respectively, for 1990-2000 (MK06).Among the inherent uncertainties in the method are those associated with fossil fuel combustion inventories and the assumed -O 2 :CO 2 stoichiometries of combustion and terrestrial photosynthesis and respiration.However, the single largest uncertainty for estimating the ocean carbon sink is the Z eff term associated with oceanic O 2 fluxes (MK06).
Early applications of the vector method assumed that the atmosphere-ocean system for O 2 was still essentially in equilibrium and thus that the net annual mean air-sea flux of O 2 was approximately zero (Keeling et al., 1993;1996).This assumption was closely related to the reasons for the decoupling of oceanic O 2 and CO 2 fluxes.First, O 2 is far less soluble than CO 2 , with only 1% of the O 2 in the oceanatmosphere system partitioning into the ocean compared to 98% for CO 2 .Second, the buffering effect of carbonate chemistry in seawater increases the air-sea CO 2 equilibration time scale by an order of magnitude relative to that of O 2 .Third, fossil fuel combustion and deforestation have raised atmospheric CO 2 by ∼35% relative to preindustrial levels, thus providing the geochemical driving force for net global oceanic CO 2 uptake.In contrast, these processes have reduced atmospheric O 2 /N 2 by only ∼240 per meg, or less than 0.03% of the total O 2 burden, since atmospheric monitoring of O 2 /N 2 began in the late 1980s (MK06).
Although the assumption of a net zero annual oceanic O 2 flux initially seemed reasonable, it has long been known that natural interannual variability in oceanic O 2 fluxes, due, e.g., to incomplete ventilation of O 2 -depleted thermocline waters in some years, followed by deep convection events in subsequent years, can create a non-zero net oceanic O 2 flux over a given year or even period of years (Keeling et al., 1993).Thus the vector method is generally applied to longer, e.g., decadal averages, of APO and CO 2 data with the assumption that the interannual variability will average to zero.In addition to natural interannual variability, recent studies have recognized that long-term ocean warming may be inducing net O 2 outgassing, with reinforcing feedbacks associated with biology and ocean stratification (Sarmiento et al., 1998;Bopp et al., 2002;Plattner et al., 2002).Both natural interannual variability and long-term net O 2 outgassing are quantified based on empirical relationships between dissolved O 2 and potential temperature extrapolated to observed heat fluxes, although these relationships are not well understood (Keeling and Garcia, 2002).MK06 estimate the combined effects of these two processes on uncertainty in ocean and land carbon sink partitioning at ±0.5 Pg C/yr.
Here we present a model-based analysis of atmospheric APO and CO 2 focusing on interannual variability in oceanic O 2 fluxes and their contribution to uncertainty in the APO vs. CO 2 method for partitioning carbon sinks.To estimate oceanic O 2 , N 2 and CO 2 fluxes, we employ a processbased ocean ecosystem model coupled to a general circulation model, which represents an advance over the OCMIPbased (i.e., phosphate-restoring) ocean carbon models that were used predominantly in previous APO studies (Stephens et al., 1998;Naegler et al., 2007).The ocean fluxes are used to force atmospheric transport model simulations that include full interannual variability in the meteorological drivers.We also force the transport model with fossil and terrestrial CO 2 fluxes, which can be scaled to O 2 based on assumed stoichiometries.Our simulation encompasses a complete set of atmospheric tracers that can be used to evaluate uncertainties in the APO vs. CO 2 method in a controlled model context in which the true surface CO 2 and O 2 fluxes are known exactly.We begin in Sect. 2 with a description of the atmospheric tracer transport model (ATM) and ocean, land, and fossil fuel fluxes as well as the APO observations used to evaluate model results.We present a comparison of model vs. observed seasonal cycles and latitudinal gradients in Sects.3.1 and 3.2, respectively.In Sect.3.3 we characterize the interannual variability in model APO and in Sect.3.4 we analyze the impact of interannual variability in oceanic O 2 fluxes on the APO vs. CO 2 method for partitioning carbon sinks.This section also includes an analysis of errors in the method associated with sparse sampling of the ATM at grid points corresponding to APO monitoring stations.We conclude with a summary of our results in Sect. 4.

Atmospheric transport model
The Model of Atmospheric Transport and Chemistry (MATCH) (Rasch et al., 1997;Mahowald et al., 1997) was used to simulate the atmospheric distribution of a range of surface O 2 , CO 2 and N 2 fluxes.The model was run at T62 horizontal resolution (about 1.9  (Kalnay et al., 1996).The surface fluxes used to force the simulations are described below and summarized in Table 1.See also Nevison et al. (2008) for additional details.
2.2 Air-sea fluxes of O 2 , CO 2 and N 2 Oceanic O 2 and CO 2 fluxes were obtained for 1979-2004 from the WHOI-NCAR-UC Irvine ocean general circulation and marine ecosystem model (abbreviated here as the WHOI model).The ecosystem model includes a nutrientphytoplankton-zooplankton-detritus (NPZD) food web with multi-nutrient limitation (N, P, Si, Fe) and specific phytoplankton functional groups (Doney et al., 1996;Moore et al., 2002Moore et al., , 2004)).The general circulation model is the NCAR ocean model, which was run at 3.6 • longitude resolution and variable latitude resolution ranging from 0.6 • near the equator to 2.8 • at mid-latitudes, and forced with historical dailyaveraged NCEP reanalysis surface wind, heat flux, atmospheric temperature and humidity, and satellite data products from 1979-2004(Doney et al., 1998, 2007)).With proper historical forcing, the NCAR ocean model captures well the magnitude and phasing of upper ocean heat content interannual variability (Doney et al., 2007).
While O 2 and CO 2 fluxes were obtained prognostically from the WHOI model, N 2 fluxes were diagnosed based on the NCEP heat fluxes (Q) as F N 2 =Q(dS N 2 /dT )/C p (Keeling et al., 1998), where dS N 2 /dT is the temperature derivative of N 2 solubility evaluated at the NCEP sea surface temperature.A similar equation was used to estimate the thermal component of the WHOI O 2 flux: F O 2 thermal =Q(dS O 2 /dT )/C p , but with two modifications recommended by Jin et al. (2007), who optimized the formula based on comparisons to explicitly modeled thermal O 2 fluxes.First, the magnitude of F O 2 thermal was scaled down by a factor of 0.7 to account for incomplete thermal equilibration of O 2 .Second, the flux was delayed for half a month to account for non-instantaneous air-sea equilibration.Similar modifications may be needed for F N 2 , but were not applied in the current study, since their appropriateness for N 2 has not yet been evaluated (X.Jin, personal communication).
In addition to the interannually-varying WHOI and NCEPbased fluxes, MATCH was also run from 1988-2004 with a variety of climatological oceanic CO 2 , O 2 , and N 2 fluxes.These included the monthly sea surface CO 2 flux climatology of Takahashi et al. (2002) (corrected to 10 m height windspeeds) and the monthly O 2 and N 2 flux anomaly climatology of Garcia and Keeling (2001) (hereafter referred to as GK01).The GK01 O 2 fluxes are based on linear regressions that use ECMWF heat flux anomalies to spatially and temporally extrapolate historical dissolved O 2 data.These empirical regressions reflect both biology and circulation and thermal-driven variations in oxygen flux.The N 2 flux anomalies are calculated from the same formula for F N 2 above, but use climatological ECMWF values rather than interannually-varying NCEP data for Q.The monthly O 2 and N 2 flux anomalies have an annual mean of zero at all gridpoints and thus do not provide realistic spatial distributions.To fill this gap, additional simulations were performed with O 2 and N 2 fluxes from the annual-mean ocean inversions of Gruber et al. (2001) and Gloor et al. (2001).The inversions are based on ocean inventory data and are subject to ocean general circulation model biases but independent of atmospheric transport influences.

Fossil fuel and land biosphere carbon fluxes
Fossil fuel carbon emissions were constructed from decadal 1 • ×1 • maps (Andres et al., 1996;Brenkert, 1998) that were interpolated between decades based on the annual regional totals of (Marland et al., 2003).The carbon fluxes were run in MATCH from 1979MATCH from -2004.The atmospheric O 2 distributions resulting from fossil fuel combustion were estimated from MATCH CO 2 multiplied by a mean fossil fuel O 2 :CO 2 stoichiometry of -1.4 (MK06).
Two versions of the CASA land biosphere model were used for terrestrial CO 2 fluxes.The first (L1 in Table 1) was the Global Fire Emissions Database (GFED v2) version of CASA, which incorporates satellite-based estimates of burned area (Randerson et al., 2005;Van der Werf et al., 2006) and for which the biomass burning component of the emissions was optimized based on CO inversion results (Kasibhatla et al., manuscript in preparation).The second (L2) was a cyclostationary (i.e., the same seasonal cycle repeated every year) "neutral biosphere" version, run in MATCH from 1979MATCH from -2004, in which annual net ecosystem production was assumed to be zero at each grid cell (Olsen and Randerson, 2004).The L1 fluxes provide more realistic interannual variability in atmospheric CO 2 than the L2 fluxes, but are limited to 1997-2004 (Nevison et al., 2008).At their original 1 • ×1 • resolution, the global mean L1 carbon flux was constrained to equal zero over the 1997-2004 period, and the L2 flux had an annual mean of zero every year.However, a small interpolation error in converting to the T62 MATCH grid resulted in a small net source of CO 2 to the atmosphere of ∼+0.1 PgC/yr for both versions of the land flux over these respective periods.

Model APO
Model atmospheric potential oxygen (APO) was calculated somewhat differently than observed APO (Eq.1), since O 2 and N 2 were run as separate tracers in MATCH.We used equation 5, in which all individual O 2 , N 2 and CO 2 tracers are in ppm of dry air and APO is in per meg units (Stephens et al., 1998;Naegler et al., 2007).
where the superscripts oc and F F denote ocean and fossil fuel tracers, respectively.X O 2 and X N 2 are the fractions of O 2 and N 2 in air (0.20946 and 0.78084, respectively).Land O 2 and CO 2 tracers could have been included in the first and third right hand side terms of Eq. 5, respectively, but simply would have canceled out, since we assumed that O land 2 = -1.1 CO land 2 .In practice, the CO F F 2 terms were important for latitudinal gradients, but had little effect on seasonal or interannual variability.For some applications, the fossil terms were omitted to yield a quantity we refer to as APO oc , which includes only the oceanic tracers.To examine the consequences of neglecting CO oc 2 in calculating APO oc , a third quantity, referred to here as O 2 /N oc 2 , was computed: For comparison to model results, we used observed seasonal cycles from Battle et al. (2006), who computed the 1996-2003 climatological cycles at 13 land-based atmospheric monitoring stations, including both the SIO (Keeling et al., 1998) and Princeton (Bender et al., 2005) networks.Battle et al. also determined the latitudinal gradient in APO from the annual mean offsets of sinusoidal fits to the observed seasonal cycles at the land-based monitoring stations and from additional ship-based measurements in the Pacific Ocean.

Seasonal cycles
The seasonal cycles in APO predicted by the WHOI-MATCH simulation (A1 in Table 1) generally agree well with the observed cycles (Figs.[1][2].The model amplitude is captured to within ±10 to 20% at most stations and the agreement in phasing is excellent at all stations, except Palmer Station, Antarctica (PSA), where the model minimum is a month  2006) observations.The angle θ from the x-axis in the polar plot is the arccosine of the correlation coefficient R between the model and observed cycle, which reflects the agreement in shape and phasing.The value on the radial axis is the ratio of standard deviations: σ model /σ obs , which represents the match between the amplitude of the model and observed seasonal cycle (Taylor, 2001).Each symbol represents 1 of the 12 stations in Fig. earlier than observed.Oceanic CO 2 , although sometimes neglected when calculating the APO seasonal cycle (e.g., GK01; Blaine, 2005), appears to contribute in a small but non-negligible manner to the amplitude agreement at some stations.Neglecting CO oc 2 leads to a greater underestimate of the seasonal amplitude at the tropical SMO and KUM stations, where the CO oc 2 and O oc 2 seasonal cycles are in phase, and a tendency to overestimate the seasonal amplitude at extratropical southern stations, where CO oc 2 opposes O oc 2 .In comparison to A1, the A3 MATCH simulation with climatological ocean fluxes yields seasonal cycles in APO that agree less well in shape and phasing with observations, and overestimate the observed amplitude by 30-40% at most extratropical stations and by nearly 100% at the Cold Bay, Alaska (CBA) station.
In general, thermal and biological O 2 fluxes tend to reinforce each other over the seasonal cycle.In fall and winter, the increased solubility of cooling surface waters causes thermal ingassing, while the breakdown of the seasonal thermocline leads to biological uptake as deep waters depleted in O 2 by microbial decomposition are ventilated (Keeling et al., 1993).Conversely, in spring and summer, warming of surface waters tends to coincide with biological O 2 outgassing due to new production.However, comparison of the WHOI model A1 and A2 simulations suggests that net total oceanic O 2 fluxes do not always correlate well with thermal O 2 fluxes (Fig. 1), as is assumed by the GK01 seasonal O 2 climatology.At southern extratropical stations in particular, thermal O 2 fluxes are only moderately well correlated in shape and phasing (R=0.4-0.9) and only account for 10-15% of the amplitude of net O 2 seasonal cycle, with the remainder due to biology and circulation.At tropical and northern stations, the thermal O 2 flux tends to correlate better and to account for a larger (20-40%) share of the seasonal amplitude.In most regions of the ocean, Figs.1-2 suggest that the explicit modeling of ocean biology and circulation provided by the WHOI model yields superior results to the climatology.Similarly, the more sophisticated ocean ecosystem dynamics of the WHOI model appears to yield more realistic seasonal cycles in APO than the P-restoring OCMIP ocean biology parameterizations tested by Naegler et al. (2007).
A caveat on the above discussion is that the choice of atmospheric transport model can influence the amplitude of the model seasonal cycle in APO (Battle et al., 2006;Naegler et al., 2007), since many atmospheric transport models produce a seasonal "rectifier" in APO (Stephens et al., 1998;Gruber et al., 2001;Blaine et al., 2005).The concept of seasonal rectification was first defined in atmospheric CO 2 studies (Pearman and Hyson, 1980;Keeling et al.,1989;Denning et al., 1995), and involves a covariation between seasonal variability in transport and surface fluxes that can yield annual-mean surface concentrations of an atmospheric tracer that are not zero, even when the annual mean surface flux of the tracer is zero.
Previous studies have found that ATMs with strong seasonal APO rectifiers, like TM3, tend to yield larger seasonal amplitudes in APO than weak rectifiers like TM2.Naegler et al. (2007) used the average of simulations with TM2 and TM3 as the best estimate of true APO, but Battle et al. (2006) suggested that TM3 may overestimate tracer concentrations and thus seasonal amplitudes due to excessive vertical trapping in surface layers.This latter suggestion is supported by the extreme overestimate of the seasonal amplitude at CBA in the A3 simulation (Figs.1-2).Transport model differences may explain the discrepancy between the results of GK01 with TM2, which captures the amplitude of the observed O 2 /N oc 2 cycle to within ±10% at most stations, and those shown in Fig. 2c, in which the A3 simulation with the same flux climatology tends to overestimate the observed APO seasonal cycle.Uncertainties associated with atmospheric transport are discussed further below in the context of latitudinal gradients.

Latitudinal gradient
In the A3 simulation, in which climatological seasonal flux anomalies have an annual mean of zero at all grid points for O 2 and N 2 , MATCH predicts non-zero annual mean APO oc .Positive values of up to ∼20 per meg are predicted over northern ocean regions centered around 45-60 • N (encompassing the CBA station) and up to ∼10 per meg at comparable southern latitudes (Fig. 3a).Similar but weaker patterns appear in the A1 simulation with seasonally varying WHOI fluxes (Fig. 3b).In contrast, the annual mean climatology simulation (A4) yields negative annual mean values of APO oc over the mid to high latitude oceans (Fig. 3c), consistent with the known patterns of net oceanic O 2 uptake in these regions.We also constructed a "composite" climatology (A5), as per Gruber et al. (2001) and Battle et al. (2006), by adding the results of MATCH simulations with the annual climatology and the GK01 seasonal anomalies.The composite climatology yields an APO oc field similar to that of the WHOI model (Fig. 3d).The results shown in Fig. 3 suggest that the MATCH transport model introduces a seasonal rectifier for simulations that include seasonally varying oceanic O 2 fluxes.Because the WHOI model itself predicts O 2 uptake over the mid to high latitude oceans (not shown), it must be that the atmospheric transport model, rather than the ocean model, is responsible for the positive APO oc values in the A1 results around 45-60 • latitude (Fig. 3b).
The MATCH seasonal rectifier in APO is caused by at least two interrelated mechanisms.First, the MATCH planetary boundary layer is shallower in summer over the ocean than in winter, due to increased convective heat release in wintertime (Fig. 4).Since positive sea-to-air O 2 flux anomalies occur during summer at mid to high latitudes, they are more likely to be trapped near the surface, resulting in a net positive APO anomaly.Stephens et al. (1998) and Gruber et al. (2001) found a similar effect with the TM2 and GCTM transport models.A second, related mechanism also contributes to the MATCH rectifier effect.The sum of convective and largescale precipitation, which is effectively a proxy for ventilation of surface layers by upward motion associated with storms, displays strong seasonal and hemispheric differences (Fig. 4).More storms occur over the ocean in the northern vs. the southern hemisphere and more storms occur in winter vs. summer in both hemispheres.Thus, negative winter O 2 fluxes are more likely to be ventilated from the boundary layer, especially in the northern hemisphere, helping to explain the winter and summer pattern in both hemispheres as well as the stronger seasonal difference in the north.
The "APO Transcom" intercomparison showed that the majority of nine transport models tested, using the GK01 seasonal anomalies to produce the quantity referred to in this study as O 2 /N oc 2 , yielded similar rectifier effects to those described above (Blaine, 2005).MATCH:NCEP fell among the strong rectifiers, but produced weaker annual mean values than some of the other ATMs in that class.TM2 fell among the weak rectifiers, especially when forced with the 1986 ECMWF winds used by GK01.
The MATCH seasonal rectifier is difficult to confirm or disprove based on comparison with observed APO.The "equatorial bulge" in APO, a prominent feature of the observations, is flattened in the A1 simulation by the strong mid to high latitude O 2 maxima associated with the seasonal rectifier (Fig. 5a,d).In contrast, the A4 simulation with the annual mean O 2 flux climatology captures the "equatorial bulge" relatively well (Fig. 5c,f).There are no seasonal rectifier effects in the A4 runs, since the O 2 fluxes are uniform throughout the year and thus can not covary with transport.The composite climatology (A5) yields a rectifier-effect-dominated latitudinal gradient that resembles the A1 results but captures the equatorial bulge better, especially when model results are subsampled in the tropical Pacific, where ship-based measurements are made, rather than zonally averaged (Fig. 5b,e).Battle et al. (2006) found similar results for TM3 simulations using the same composite climatology.
Aside from the equatorial bulge, both the A1 and A5 simulations produce a south-to-north gradient that agrees reasonably well with observations (Fig. 5d,e).Both simulations, particularly A5, overestimate the observed annual mean value at CBA, suggesting an exaggerated seasonal rectifier effect in that area of the North Pacific.However, the annual mean simulation (A4), which lacks a seasonal rectifier, substantially underestimates APO at midlatitudes relative to the tropics, particularly in the northern hemisphere.Collectively, the above results suggest, in agreement with Naegler et al. (2007), that the APO seasonal rectifier probably exists to some degree in the real world, although both MATCH and TM3 may tend to exaggerate it.Overall, our results support the conclusion of Naegler et al. (2007) that evaluating ocean models (or flux climatologies) based on the latitudinal gradient in APO is compromised by uncertainties in atmospheric transport models.The interannual variability and complexity of the observations further complicate the evaluation (Battle et al., 2006).

Interannual variability
A breakdown of APO from the WHOI-MATCH simulation (A1) into its components shows that oceanic O 2 fluxes are the main drivers of interannual variability in model APO (Fig. 6a,b).Most of the variability originates at mid to high latitudes, especially in the Southern Ocean (Fig. 6a).CO oc  the relatively low total variability predicted there (Fig. 6b).
In general, interannual variability in the MATCH simulation with WHOI total O 2 (W2 in Table 1) is moderately correlated to variability in thermal O 2 (W4) (R∼0.6 to 0.8), but the ratio of total:thermal variability is spatially heterogeneous, with highest ratios in the Southern Ocean and lowest ratios in the tropics (Fig. 6c).The underlying total O 2 flux:heat flux ratios of the WHOI model, which drive the patterns in Fig. 6c, range from >10 nmol/J in some seasons in the Southern Ocean to 1-2 nmol/J in the tropics and <0 in the equatorial upwelling zone, where, in contrast to most ocean regions, solubility and biology and circulation act in opposite directions on O 2 (Stephens et al., 1998).These patterns are generally consistent with O * 2 vs. θ slopes observed in the ocean thermocline (Keeling and Garcia, 2002), where O * 2 is O 2 corrected for biological activity based on the observed phosphate concentration and the Redfield O 2 :P ratio (e.g., Gruber et al., 2001).
The amplitude of interannual variability in simulated APO ranges from ±5 to ±10 per meg year −1 at the 12 monitoring stations in Fig. 8. Observed APO from the Princeton University network (Bender et al., 2005) shows a similar to somewhat larger range of variability (Fig. 8).APO observations from the SIO network (not shown in Fig. 8) vary on the order of 5 to 10 per meg year −1 (Hamme and Keeling, 2008 1 ).The A1 model results show higher frequency variability than the Princeton data (Fig. 8), probably due, at least in part, to the different smoothing algorithms applied.The Princeton data have been smoothed and deseasonalized with a 620-day smoothing algorithm (Thoning et al., 1989), whereas we have used a 13-month weighted running average for the model output.The SIO data, which are deseasonalized with a four-harmonic fit and smoothed with a six-month Pole and expressed in per meg units.Heavy black line is total APO, including fossil fuel; solid gray line is oceanic CO 2 ; dotted black line is oceanic O 2 ; dashed black line is oceanic N 2 ; dash-dot black line is fossil fuel.running mean, also show higher frequency variability than the Princeton data (Hamme and Keeling, 2008 1 ).In general, the A1 simulation appears to predict a realistic range of variability in APO (±5 to ±10 per meg year −1 ), although model APO is not well correlated temporally with either the Princeton (Fig. 8) or SIO observations (Hamme and Keeling,  2008 1 ).

Partitioning carbon sinks using APO and CO 2
The MATCH simulations described in Table 1 offer a selfcontained system for exploring some of the uncertainties in the APO-vs.-CO 2 method for partitioning land and ocean carbon sinks in a more controlled manner than is possible in the real world.The model simulation contains no measurement error or drift in CO 2 / t and APO/ t and no uncertainties in α bio , α f , F fuel , or Z eff .In addition, F land and F ocean can be computed exactly by globally integrating www.biogeosciences.net/5/875/2008/Biogeosciences, 5, 875-889, 2008  the fluxes used to force MATCH.These "true" fluxes can be compared to the F land and F ocean terms inferred by solving Eqs. 2 and 4 using the time derivatives of MATCH APO and total CO 2 (Fig. 9).Having eliminated most conventional sources of uncertainty, the model analysis allows us to examine the uncertainty in the method arising from 1) variability introduced by transport acting on surface fluxes and sampling of the resulting atmospheric tracers at a limited number of stations, and 2) the O 2 outgassing term Z eff , which, like the true carbon sinks, can be calculated exactly by globally integrating the oceanic O 2 fluxes used to force MATCH.Early applications of the vector method solved the O 2 /N 2 mass balance (Eq.3) first for F land and then used that result to solve the CO 2 mass balance (Eq.2) for F ocean (Keeling et al., 1996).We have followed the more recent preference for solving the APO mass balance (Eq.4) first for F ocean and then using Eq. 2 for F land (Bender et al., 2005;MK06).This second approach offers several advantages over the O 2 /N 2 vs. CO 2 system, including the ability to use "global" atmospheric CO 2 datasets from the NOAA CCGG network to solve Eq. 2 instead of CO 2 data from the more limited number of O 2 /N 2 monitoring stations (MK06).We calculate APO/ t from the A1 simulation and CO 2 / t by adding either W1+FF+L1 or W1+FF+L2.
Although our ultimate goal is to examine the uncertainty in the method associated with O 2 outgassing, we first examine the method's performance when the exact Z eff is used in Eqs. 2 and 4. In these first tests, any discrepancies between the estimated and true F land and F ocean terms must be introduced by atmospheric transport and sampling errors.Figure 10a demonstrates that transport and sparse spatial sampling introduce large uncertainties, especially in F land , when individual station data are used to estimate CO 2 / t.This transport and sampling uncertainty is reduced when global CO 2 / t (calculated by integrating all gridpoints in the MATCH surface layer) is substituted for the individual station data (Fig. 10b).The uncertainty is further reduced when APO/ t is calculated from the average of a group of stations rather than a single station (Fig. 10c).In this figure, calculations with L1 (with realistic interannual variability in land CO 2 ) yield larger deviations from true F land than calculations with L2, the cyclostationary land CO 2 tracer.All of the station groupings shown in Fig. 10c, which were chosen because they have been used in real-life applications, give a good approximation of the truth.The ALT, LJO and CGO combination (favored by MK06) is used in Fig. 10d to explore the impact of staggered decadal intervals on the APO vs. CO 2 method.Figure 10d shows that staggered decadal time spans (using L2) have some scatter associated with transport and sampling uncertainty, but give true F land and F ocean to within ±0.2 PgC/yr or better, with a standard deviation of only ∼0.05 PgC/yr.From the above results, we conclude that transport and sparse sampling introduce a small error of ∼±0.2 PgC/yr into the CO 2 land and ocean sinks inferred using the APO vs. CO 2 method.This error reflects uncertainty not only in the partitioning between land and ocean sinks, but also where f O oc 2 and f N oc 2 are the globally integrated WHOI ocean fluxes in moles/yr used to force MATCH, M=1.768×10 20 is the total moles of dry air in the atmosphere, and β, γ and α bio are conversion factors defined in the text that convert per meg to Pg C. Gray italic text denotes quantities known in the model simulation but unknown in the real world.Unless labeled otherwise, all these quantities are in PgC/yr, averaged over the 8 year interval 1997-2004.Note that F land in the model is actually a small source of carbon and thus has a negative value in the diagram, since F land and F ocean are defined as positive when they act as sinks.
in the absolute magnitude of the combined sinks (which is set by the difference F fuel − CO 2 / t), since the transport and sampling uncertainty affects CO 2 / t.The analysis in Fig. 10 does not necessarily identify an optimal spatial sampling strategy, since the results may be influenced by transport biases in the particular atmospheric transport model (MATCH) used here.It is possible, for example, that spurious seasonal rectifier effects for both APO, as discussed above, and CO 2 (Denning et al., 1995;Gurney et al., 2003) may bias the model uncertainty relative to the real world.
In the results discussed thus far, the true Z eff is known exactly and used to solve Eqs. 4 and 2. Figures 10c,11 and 12 investigate what happens in the more realistic case where the O 2 outgassing term is not exactly known.Figure 10c shows that ignoring O 2 outgassing entirely leads to a clear underestimate of the ocean carbon sink F ocean .This is true whether the L1 or L2 land CO 2 tracer is used.The neglect of Z eff leads to a compensating overestimate of the true land sink in the L2 case, but partly cancels transport and sampling error in F land in the L1 case, such that the true F land is estimated about equally well with or without Z eff (Fig. 10c).
In real-life applications of the APO vs. CO 2 method, O 2 outgassing is not ignored (Battle et al., 2000;Bender et al., 2005;MK06).Rather, current studies have an arguably good estimate of the long-term mean Z eff (Bopp et al., 2002;  Keeling and Garcia, 2002;Plattner et al., 2002).However, these studies still assume that interannual variability in oceanic O 2 fluxes averages to zero for the time period (commonly ∼1decade) over which the method is applied.Figure 11 shows that the decadal assumption is reasonable for the WHOI model.In Fig. 11a, four different 10-year average O 2 fluxes, calculated at staggered intervals over the 25-year span of the W2 simulation, deviate relatively little (σ =±6 Tmol O 2 /yr) from the long-term 25-year mean of 30 Tmol O 2 /yr.This is true despite the relatively large range of interannual variability (-50 to +85 Tmol O 2 /yr) in the O 2 fluxes.The model 25-year mean corresponds closely to the value of Z eff used in calculations involving observed APO and CO 2 (Keeling and Garcia, 2002;MK06).When this 25year mean is substituted for the true 10-year averages in the APO vs. CO 2 mass balance calculations, the resulting F ocean deviates only slightly more from the truth than the F ocean estimated when the exact Z eff is used, given the small transport and sampling-related background error that already exists in the calculation (Fig. 11b).
When one moves to shorter, e.g., 5-year, averages, the above conclusions start to break down.The eight staggered 5-year average O 2 fluxes shown in Fig. 11a deviate more substantially (σ =±16 Tmol O 2 /yr) from the 30 Tmol/yr longterm mean than did the 10-year averages.Accordingly, the values of F ocean calculated by substituting the long-term mean Z eff for the true 5-year averages display larger deviations from the truth, with errors ranging up to ±0.5 Pg C/yr.For even shorter time intervals, e.g., 1-year averages, the problem is compounded further, leading to large errors and unrealistic year-to-year swings in the inferred value of F ocean of ±1 PgC/yr or more, consistent with previous studies (McKinley et al., 2003;Bender et al., 2005).
To quantify more systematically the error in the APO vs. CO 2 method as a function of time span, we solve Eqs. 4 and 2 using the MATCH A1, and W1+FF+L2 time series over t intervals ranging from 1 to 15 years and compute the standard deviation from the "truth" in F ocean .We consider all whole-year staggered intervals for each respective t between 1980 and 2003.We perform these calculations using both the long-term mean Z eff and the true Z eff in order to distinguish between errors due to uncertainty in O 2 outgassing versus in transport and sampling.For calculations using true Z eff , the standard deviation σ of F ocean from the truth is <0.1 PgC/yr for t=3 to 15 years, with a small increasing trend in σ as t decreases (Fig. 12).Intervals of t≤3 years show a steeper increase in σ associated with the increased uncertainty in APO/ t and, to a lesser extent, global (surface level) CO 2 / t.This result suggests that time spans of 3 years or less are too short to yield accurate estimates of F ocean due to errors associated with sparse sampling and transport, even if O 2 outgassing is known exactly.Calculations substituting the 1980-2003 mean Z eff for the true Z eff yield a somewhat larger deviation from the true F ocean for all time intervals, but σ remains <0.1 PgC/yr for t ranging from 15 to ∼9 yrs.However, as t decreases to 9 years or less, σ begins to increase more sharply, paralleling the increasing standard deviation of true Z eff from long-term mean Z eff (Fig. 12).For t≤3 years, the increase in σ steepens still more, due to the increased influence of transport and sampling uncertainties, as discussed above.

Conclusions
Forward simulations with the MATCH atmospheric transport model, using surface fluxes from the WHOI ocean ecosystem model, appear comparable or superior to those with observed climatological ocean fluxes in terms of reproducing observed seasonal and spatial patterns in atmospheric potential oxygen (APO).However, atmospheric transport uncertainties, especially those associated with the seasonal rectifier effect in APO, compromise our ability to evaluate ocean models Fig. 12. Summary of standard deviations of inferred minus true F ocean (in Pg C/yr).Calculations use deseasonalized WHOI-MATCH results over t intervals ranging from 1 to 15 years, in which all whole-year staggered intervals for each respective t between 1980 and 2003 are considered (e.g., N=24 intervals for t=1 and N=10 intervals for t=15).Black solid line shows results using true Z eff over the interval t; black dotted line shows results using the 1980-2003 mean Z eff ; Gray dot-dash line shows the standard deviation of mean minus true Z eff (converted to PgC/yr as described in Fig. 9) for the N staggered intervals at each t.All calculations use global CO 2 / t from the W1+FF+L2 simulation and the ALT, LJO, CGO station average for APO/ t from the A1 simulation.
with APO.Oceanic O 2 fluxes dominate seasonal, spatial and interannual variability in simulated APO, but oceanic CO 2 fluxes can contribute substantially to all three types of variability, especially in the tropics.The WHOI model predicts strong interannual variability in the annual-mean ocean O 2 flux (range -50 to +85 Tmol O 2 /yr), which is concentrated largely at high latitudes.An analysis of the APO vs. CO 2 mass balance method for partitioning land and ocean carbon sinks, performed in the controlled context of the MATCH simulation, in which the true surface carbon and oxygen fluxes are known exactly, suggests a small uncertainty (ranging up to ±0.2 Pg C) in the sinks due to errors associated with atmospheric transport and sparse sampling.Natural interannual variability in ocean O 2 fluxes introduces only a small additional error into the inferred carbon sink partitioning when the APO-vs.-CO 2 method is applied over decadal periods, but this error becomes larger when the method is applied over shorter intervals.

Fig. 2 .
Fig. 2.Taylor diagrams summarizing the comparison of MATCH seasonal cycles toBattle et al. (2006) observations.The angle θ from the x-axis in the polar plot is the arccosine of the correlation coefficient R between the model and observed cycle, which reflects the agreement in shape and phasing.The value on the radial axis is the ratio of standard deviations: σ model /σ obs , which represents the match between the amplitude of the model and observed seasonal cycle(Taylor, 2001).Each symbol represents 1 of the 12 stations in Fig.1, or a 13th station, Amsterdam Island, France (AMS).The 13 stations are sorted into 4 latitude bands according to the figure legend.The following MATCH simulations are shown: (a) WHOI APO oc (A1) (b) WHOI O 2 /N oc 2 (A1, but CO oc 2 not included), (c) seasonal climatology APO oc (A3).

Fig. 4 .
Fig. 4. Zonal averages (over ocean regions only) exploring the cause of the seasonal rectifier effect for (a) January, (b) July.Black solid line is the WHOI oceanic O 2 flux (mol m −2 yr −1 ); black dotdash line is the WHOI-MATCH O 2 seasonal anomaly (ppm) from the W2 simulation; gray dotted line is the MATCH planetary boundary layer height (pblh); gray solid line is the sum of MATCH large scale and convective precipitation (prect), a proxy for ventilation of surface layers by upward motion associated with storms.Both pblh and prect are normalized at each latitude band according to the annual mean value in that band.Prect is masked equatorward of 30 • because convective precipitation has strong patterns there that distract from the mid-latitude patterns the figure seeks to emphasize.

Fig. 5 .
Fig. 5. Left column: Components of MATCH zonally averaged annual mean APO, in which all values are normalized to the South Pole and expressed in per meg units.Heavy black line is total APO, including fossil fuel; solid gray line is oceanic CO 2 ; dotted black line is oceanic O 2 ; dashed black line is oceanic N 2 ; dash-dot black line is fossil fuel.(a) WHOI (A1) (b) composite climatology (A5), (c) annual climatology (A4).Right column: Compares MATCH annual mean latitudinal gradient in APO to Battle et al. (2006) station (gray squares) and shipboard (gray triangles) data.Black squares are MATCH annual mean values at the Battle et al. (2006) stations, which differ in some cases from the zonal averages (thin black line).Black triangles are MATCH annual-mean values along the Pacific transect described by Battle et al. (2006).Since the values at the SPO are arbitrary, all results have been shifted to aid visual comparison.
Fig. 5. Left column: Components of MATCH zonally averaged annual mean APO, in which all values are normalized to the South Pole and expressed in per meg units.Heavy black line is total APO, including fossil fuel; solid gray line is oceanic CO 2 ; dotted black line is oceanic O 2 ; dashed black line is oceanic N 2 ; dash-dot black line is fossil fuel.(a) WHOI (A1) (b) composite climatology (A5), (c) annual climatology (A4).Right column: Compares MATCH annual mean latitudinal gradient in APO to Battle et al. (2006) station (gray squares) and shipboard (gray triangles) data.Black squares are MATCH annual mean values at the Battle et al. (2006) stations, which differ in some cases from the zonal averages (thin black line).Black triangles are MATCH annual-mean values along the Pacific transect described by Battle et al. (2006).Since the values at the SPO are arbitrary, all results have been shifted to aid visual comparison.

Fig. 6 .
Fig. 6.(a) Interannual RMS variability 1988-2004 in WHOI-MATCH APO oc (A1 simulation) computed as the RMS of the differences between each month and the corresponding month from the MATCH climatological seasonal cycle (e.g., model January 1998 minus model climatological January).(b) Ratio of RMS variabilities 1.1 CO oc 2 :O oc 2 (W1:W2).(c) Ratio of RMS variabilities O oc 2 :O thermal 2

Fig. 7 .
Fig. 7. Spatial coherence of interannual variability in WHOI-MATCH APO (A1 simulation).Shows the correlation coefficient R between the Cape Grim, Tasmania station and all other gridpoints in the 1988-2004 deseasonalized A1 time series.Station location is indicated by a blue X.

Fig. 8 .
Fig. 8. Interannual variability in APO from the A1 WHOI-MATCH simulation at the stations (a) to (l) defined in Fig. 1.Model APO is calculated from monthly mean output using Eq. 5.The seasonal cycle is removed with a 13-month weighted running average and interannual variability in per meg year −1 is calculated as the time derivative of the deseasonalized time series.Observed interannual variability from Bender et al. (2005) is shown (blue curve) at BRW, CGO, SMO, and SYO.

Fig. 9 .
Fig. 9. Vector diagram showing MATCH time derivatives APO/ t (from A1 simulation) and CO 2 / t (from W1+FF+L1 simulations) sampled at station Alert.The carbon sinks F land and F ocean inferred by solving Eqs. 4 and 2 are compared to the "actual" surface fluxes used to force the MATCH simulation.The time derivative of the MATCH fossil fuel tracer CO F F 2 / t (simulation F F ) is also compared to the actual surface flux F fuel used to force MATCH.The O 2 outgassing term (Z eff ) is computed as

Fig. 10 .
Fig. 10.Summary of calculations testing the sensitivity of F land and F ocean , estimated by solving Eqs. 4 and 2 using WHOI-MATCH tracers, to different spatial sampling strategies, time spans, land CO 2 tracers, and inclusion or omission of the O 2 outgassing term Z eff .APO/ t is calculated from the A1 WHOI-MATCH simulation and CO 2 / t is calculated by adding either W1+FF+L1 or W1+FF+L2.The difference between inferred and actual F land and F ocean is shown on the x and y-axis, respectively, in Pg C/yr.The "truth" at x=0, y=0 is marked with a star.(a) 1997-2004 using the L1 land CO 2 tracer.APO/ t and CO 2 / t are evaluated separately at 13 individual stations, as labeled on plot, using legend symbols described in Fig. 2. (b) same as (a) except using global CO 2 / t.(c) 1997-2004 using L1 (circles) or L2 (squares) land CO 2 tracer.The numbers 1-5 indicate the combination of stations used to calculate average APO/ t: 1=ALT, LJO; 2=ALT, LJO, CGO; 3=BRW, SMO, CGO; 4=ALT, LJO, KUM, CGO, PSA; 5=ALT, CBA, LJO, KUM, SMO, CGO, PSA.Global CO 2 / t is used for all cases.Gray symbols are corresponding solutions when the O 2 outgassing correction is omitted.(d) Shows 17 staggered 10-year intervals, using the L2 land CO 2 tracer, from 1979-1988 (lightest circle) to 1995-2004 (darkest).Uses global CO 2 / t and the ALT, LJO, CGO station average for APO/ t.The 4 intervals featured in Fig. 11 are explicitly labeled.
• E) vs. MQA (54 • S, 159 • E) is ∼0.65.Similar correlation coefficients are found among the SIO networks stations within each hemisphere, suggesting the presence of common largescale signals as well as local variability around each station