Regional air-sea fluxes of anthropogenic carbon inferred with an Ensemble Kalman Filter

[ 1 ] Regional air-sea fluxes, ocean transport, and storage of anthropogenic carbon ( C anth ) are quantified. Observation-based C anth data from the ocean interior are assimilated into the Bern3D dynamic ocean model using an Ensemble Kalman Filter. Global uptake of C anth is estimated to be 131 ± 18 GtC over the period 1770 to 2000. Uncertainties from systematic biases in the reconstruction of C anth are assessed by assimilating data from four global and six Atlantic reconstructions and found to be comparable or larger than uncertainties from ocean transport. Aggregated fluxes for the southern high-latitude, tropical and midlatitude, and northern high-latitude ocean agree within 0.11 GtC a (cid:1) 1 for the two reconstructions with the highest skill score, whereas regional uptake rates are up to a factor of three different. Results indicate that uptake and regional partitioning of anthropogenic carbon in the Southern Ocean remains uncertain.


Introduction
[2] A quantitative understanding of the uptake of anthropogenic carbon, C anth , by the ocean is required to establish the atmospheric budget of carbon dioxide (CO 2 ) and to project the future evolution of CO 2 and climate. The air-sea fluxes of carbon are currently perturbed by rising atmospheric CO 2 resulting from carbon emissions from fossil fuel burning and land use. Unfortunately, it is not possible to measure the ocean uptake flux of anthropogenic carbon directly, as the unperturbed (preindustrial) pattern of the airsea exchange is not well known. Therefore, there are large uncertainties in the spatial pattern and magnitude of the airsea flux of C anth . The purpose of this study is to quantify regional air-sea fluxes of C anth and their uncertainties. We apply an Ensemble Kalman Filter (EnKF) [Evensen, 2003] to assimilate data of the anthropogenic carbon concentration from the ocean interior into the Bern3D model [Müller et al., 2006], thereby constraining air-sea fluxes at the surface and meridional transport rates at depth.
[3] The distribution of anthropogenic carbon has to be reconstructed from measurements of total inorganic carbon (C T ) and other tracers [Chen and Millero, 1979;Gruber et al., 1996;Touratier and Goyet, 2004;Touratier et al., 2007;Waugh et al., 2004;M. Vázquez-Rodríguez et al., Reconstructing preformed properties and air-sea CO 2 disequilibria for water masses in the Atlantic from sub-surface data: An application in anthropogenic carbon determination, submitted to Journal of Marine Systems, 2008, under review]. This is usually done by estimating the preindustrial, steady-state distribution of C T using measurements of a range of tracers in waters that are known to be unaffected by the anthropogenic CO 2 perturbation. A main goal of this work is to quantify uncertainties arising from the reconstructions of anthropogenic carbon in the ocean. Data from six different reconstruction methods are assimilated into the model. These methods are the DC* method [Gruber et al., 1996], the CFC-shortcut (or CFC) method [Thomas and Ittekkot, 2001], the TTD method [Waugh et al., 2006], the TrOCA method [Touratier and Goyet, 2004;Touratier et al., 2007], the 8 À C T 0 method (M. Vázquez-Rodríguez et al., manuscript under review, submitted 2008) and the Institut Pierre Simon Laplace (IPSL) method [Lo Monaco et al., 2005]. A short introduction to these different C anth reconstruction methods is given in section 2. Different inverse methods have been employed in carbon cycle research to quantify variables that cannot be measured easily, such as detecting sources and sinks of atmospheric CO 2 [Enting and Mansbridge, 1989;Tans et al., 1990], export fluxes in the ocean with an adjoint model [Schlitzer, 2002], air-sea exchange of anthropogenic carbon [Gloor et al., 2003;Mikaloff Fletcher et al., 2006] or estimating different surface trace gas fluxes [Heimann and Kaminski, 1999].
[4] The Ensemble Kalman Filter is a sequential data assimilation technique where observations are assimilated into the model whenever they are available. The EnKF has gained popularity in recent years because of its simple conceptual formulation and relative ease of implementation [Evensen, 2003], as well as applicability to many different problems. In the EnKF, an ensemble of model simulations is performed and the error covariance of the model state is minimized compared to observations. Typical EnKF appli-cations are the assimilation of spatiotemporally varying data into dynamical models to constrain model states in prediction systems such as those used for weather forecast. The EnKF can also be applied to equilibrium model solutions; for example, Hargreaves et al. [2004] determined parameters of a dynamical ocean model by assimilating the distribution of temperature and salinity in an EnKF scheme. Recently, Mikaloff Fletcher et al. [2006] assimilated ocean interior data of anthropogenic carbon, C anth , to determine regional air-sea fluxes and meridional transport rates of anthropogenic carbon. Mikaloff Fletcher et al. [2006] relied on C anth data from GLODAP [Key et al., 2004], where C anth was estimated from measurements of dissolved inorganic carbon and other biogeochemical tracers using the DC* reconstruction method developed by Gruber et al. [1996]. Mikaloff Fletcher et al. [2006] applied pulse response functions from a suite of models, including the Bern3D model, which is used in this study, to compute rates of surface-to-deep transport of C anth and optimized airsea fluxes by minimizing deviations between modeled and reconstructed C anth . Errors in data-based estimates of C anth have been addressed by investigating possible biases in the DC* reconstruction [Matsumoto and Gruber, 2005]. Here, we extend this earlier work (1) by assimilating C anth data from six different, published reconstruction methods to quantify uncertainties related to the C anth input data and (2) by applying the EnKF data assimilation approach to evaluate the sensitivity of the results to the inversion method and to test the general applicability of the EnKF-Bern3D method for biogeochemical applications.
[5] The identification and quantification of uncertainties is important for any inverse analysis. In the ocean case, errors in model transport, which are due to limited model resolution and imperfect parametrization of physical and biochemical processes, appear to be one of the biggest sources of uncertainty [Gruber et al., 2001;Mikaloff Fletcher et al., 2006]. Another source of uncertainty is the uncertainty of the databased estimates that constrain the optimization. Measurements errors related to precision are normally included in optimizations, but these random errors do not account for systematic biases in the data. However, systematic errors may affect results considerably, and thus should be quantified.
[6] We address these two sources of uncertainty with emphasis on uncertainties from C anth data. This allows us to compare the magnitude of the two major sources of uncertainty and to quantify how strong biases in ocean model transport and in estimates of C anth affect the inferred air-sea fluxes and meridional transport rates. We note that Mikaloff Fletcher et al. [2006] thoroughly assessed the influence of different model biases by applying different ocean models. Here, the Bern3D ocean model is applied with different circulation patterns to test whether resulting uncertainties are comparable and sensitivities similar as for the suite of models applied by Mikaloff Fletcher et al. [2006]. The second source of uncertainty is addressed by assimilating data-based estimates from six different anthropogenic CO 2 reconstructions methods for the Atlantic and from four different methods for the global ocean. Thus, uncertainties related C anth data are assessed in a much more thorough way than it has been done previously. The air-sea fluxes in the Atlantic are constrained with C anth data from four sections calculated with all of the six reconstruction methods. For the global ocean, gridded data obtained with the DC*, CFC, TTD and TrOCA methods are used.

Data
[7] For global analysis and comparison with the study of Mikaloff Fletcher et al. [2006], section data from the GLODAP database [Key et al., 2004] are used. C anth based on the DC* [Gruber et al., 1996], the CFC [Thomas and Ittekkot, 2001] and the TTD [Waugh et al., 2004] methods are provided on the GLODAP site. In addition, we estimated C anth with the TrOCA method [Touratier et al., 2007] from the GLODAP section data using carbon, oxygen, alkalinity (Alk) and temperature. Additionally, data from four different sections and from six different reconstruction methods [Vázquez-Rodríguez et al., 2008] are applied to infer fluxes in the Atlantic and to further explore and quantify the impact of systematic uncertainties in input data. Estimates from the 8 À C T 0 (M. Vázquez-Rodríguez et al., manuscript under review, submitted 2008) and the IPSL [Lo Monaco et al., 2005] method are available for these sections in addition to the four methods used in the global analysis. The sections are WOCE A14 (1995), WOCE I06-Sb (1996), CLIVAR Repeat Section A16N (Legs 1, 2) (2003) and WOCE AR01 (1998) (Figure 1). These four sections include the Atlantic water masses from 60°S to 60°N and are located in or close to (WOCE I06-Sb) the Atlantic basin.
[8] All available data points from the different sections are individually assigned to the appropriate grid cell of the Bern3D model. In the EnKF assimilation, the ocean model results are sampled for each observation-based data point at the year when the water sample was obtained and from the grid cell associated with the data point. This ensures that we fully account for the spatiotemporal evolution of the input data. The methods used to reconstruct C anth are briefly described here; the reader is referred to the original publications for further details. We can classify the six different reconstruction methods in two groups; the first (DC*, TrOCA, 8 À C T 0 and IPSL) includes back calculation techniques, which allows one to separate C anth from the measured C T in the water column, including C T variations due to circulation, remineralization and dissolution of calcium carbonate. The second group (CFC and TTD) contains transient tracer-based methods which rely on estimating the age of a water parcel. These methods estimate C anth without direct measurements of carbon.
[9] The DC* method [Gruber et al., 1996] is based on the quasi conservative tracer C* and the assumption that water masses exchange along isopycnal surfaces. C* accounts for the contributions of calcium carbonate dissolution and remineralization of organic matter to the measured dissolved inorganic carbon concentration (C) by assuming that these fluxes are linked to oxygen (O 2 ) and Alk through constant and uniform Redfield ratios (r) (C* = C À r C:O2 O 2 À 1 2 (Alk + r N:O2 O 2 )). The DC* term is defined as the difference between present C* and preformed C* ,0 and represents the sum of C anth and a contribution from pCO 2 disequilibrium between air and sea (DC*(t) = C*(t) À C* 0 = C anth + DC dis ). Preformed C* ,0 is assumed to be time-invariant. Preformed oxygen, O 2 0 , required to compute preformed C*, is estimated by assuming that O 2 0 equals the concentration in equilibrium with the atmosphere at in situ temperature and salinity. Preformed alkalinity, Alk 0 , is computed from in situ tracer data by applying a multiple linear regression parameterization. The parameterization has been derived from surface water data. Finally, the DC* method considers a constant CO 2 air-sea disequilibrium, DC dis that is subtracted from the DC* tracer to get an estimate of C anth . DC dis is estimated for each isopycnal layer either from old C anth free waters for deep-reaching isopycnals or alternatively by assigning a single water mass age to each water parcel using CFC or other tracers with a time information.
[10] The CFC-shortcut method [Thomas and Ittekkot, 2001], only uses CFC concentration measurements. A single time, t, denoting when a water parcel (sample) has lost contact with the atmosphere, is calculated from the measured CFC concentration. Carbonate chemistry equations are then used to compute the C T concentration in equilibrium with the known atmospheric CO 2 partial pressure at time t and for the preindustrial time period. C anth is then the difference between these two equilibrium C T concentrations. In contrast to the approaches above, no assumptions about the Redfield ratios between the different biogeochemical elements are required. It is implicitly assumed that a water parcel is characterized by a single discrete age, that the air-sea disequilibrium has not changed since preindustrial time, and that ocean circulation has remained in steady state. The Transit Time Distribution (TTD) method [Waugh et al., 2006[Waugh et al., , 2004 can be viewed as a sophistication of the CFC-shortcut method. A distribution for the transit time is assigned to a sample rather than a discrete age as in the CFC-shortcut method, thereby accounting for the mixing of water masses with different histories. The time distribution is based on the measurements of different transient tracers such as CFCs or tritium and helium.
[11] The TrOCA method [Touratier and Goyet, 2004;Touratier et al., 2007], uses a quasi-conservative tracer, similar to those of NO and PO [Broecker, 1974], called TrOCA which is based on the Redfield equation. In order to determine the anthropogenic carbon concentration in sea- water, the conservative tracer TrOCA 0 (TrOCA prior to anthropogenic contamination), is subtracted from the quasiconservative tracer TrOCA. Since TrOCA 0 can simply be calculated using the measured temperature and total alkalinity, anthropogenic carbon concentrations can be directly determined with measurements of oxygen, total alkalinity, total inorganic carbon and temperature. The 8 À C T 0 method (M. Vázquez-Rodríguez et al., manuscript under review, submitted 2008) constitutes a revised version of the C* method. The main modifications are (1) DC dis and preformed total alkalinity A T 0 are parameterized taking subsurface layer (100 -200 meter depth) data as the only reference and assuming a multi-end-member mixing model instead of an isopycnal mixing model; (2) the air-sea CO 2 disequilibrium DC dis is not taken as constant over time, but its spatial variability and temporal evolution since preindustrial time is modeled; (3) no age estimates from CFCs or other age tracers are necessary for calculating anthropogenic carbon, except for the case of subsurface reference samples.
[12] In the IPSL method [Lo Monaco et al., 2005], anthropogenic carbon is estimated using the preformed C T method, a back-calculation technique developed by Brewer [1978] and Chen and Millero [1979]. The method has been tested many times, e.g., in the North Atlantic [Körtzinger et al., 1998] and in the Southern Ocean [Lo Monaco et al., 2005]. In the latter study, the method was modified to allow oxygen disequilibria with the atmosphere on the basis of the DC* method by Gruber et al. [1996], and an optimum multiparameter analysis was used to determine the relative contributions of northern and southern waters to the observed properties. Substantial differences in the reconstructed column inventories are found among methods [Vázquez-Rodríguez et al., 2008] (Figure 1). For example, the reconstructed C anth is only 14 mol m À2 at 60°S for the DC* method, but around 50 mol m À2 for the 8 À C T 0, the TrOCA and TTD method, and around 75 mol m À2 for the IPSL method, and even higher for the CFC-shortcut method. Such large differences of a factor of three are quite astonishing, as computational differences between the reconstruction methods are often quite subtle.
[13] In general, the CFC-shortcut method yields much higher column inventories than all of the other methods ( Figure 1). We attribute this to the application of a single time to characterize the age of a water parcel. This assumption does not seem very realistic and results from this method should be viewed with some caution [Matear et al., 2003;Waugh et al., 2006]. The IPSL approach yields somewhat higher column inventories along the four sections than the remaining four methods, except in high northern latitudes of section A16N. Consequently, we expect that the assimilation of the IPSL data will yield a high estimate for the air-sea flux of anthropogenic carbon in the Atlantic. The TrOCA and 8 À C T 0 methods yield very similar column inventories for all four sections, while inferred inventories from the TTD method are smallest along A16N (11°S to 66°N), but usually larger than the TrOCA and 8 À C T 0 values for I06-Sb (66°S to 28°S).
[14] There are also differences in the magnitude and sign of the meridional gradients in C anth (Figure 1). Differences are particularly large for the Southern Ocean. For example, the DC* data suggest increasing column inventories from 66°S to 40°S (section I06-Sb), whereas the IPSL data suggest a decreasing trend. In conclusion, the different reconstruction methods yield substantially different column inventories and meridional gradients in C anth when applied to the same section data. Consequently, we expect different solutions for the air-sea flux and meridional transport of C anth when assimilating C anth data from the different reconstruction methods.

Model Setup 2.2.1. Ensemble Kalman Filter
[15] The EnKF is a sequential filter method to estimate the ''true'' state of a model, by assimilating observations into the model. In general application, the model ensemble is integrated forward in time, and whenever measurements are available, these are used to update the model states before the integration is continued. The update is performed in a way that minimizes the root mean square difference between the data and the model forecast.
[16] A brief explanation of how we implement the EnKF method is given below; for a more thorough description of the EnKF see Evensen [2003Evensen [ , 2004. A model forecast, i.e., an estimated model state, y m f , is introduced for each ensemble member m 2 1, . . ., N. Here, a set of 22 air-sea flux scaling parameters, y m f (l), (l 2 1, . . ., L = 22) is prescribed for 22 different regions in the Bern3D model. The ensemble includes N = 32 members, which are iteratively optimized. An iterative step includes (1) an ensemble of 32 simulations over the industrial period with the Bern3D model, each run with a different set of 22 prescribed air-sea gas exchange scaling parameters, y m f ; and (2) optimized or analyzed flux scaling parameter, y m a , are computed with the EnKF optimization scheme.
[17] After each optimization step, the models are reinitialized with the corresponding analyzed parameter y m a ; in other words, the forecast for the next iteration is set equal to the analyzed states: y m f = y m a of the previous iteration. Then, step 1 and 2 are repeated until convergence of the solution.
[18] The ensemble of analyzed model states y a = (y 1 a , . . ., y N a ) is computed from the ensemble of forecast, y f = (y 1 f , . . ., y N f ), the assimilated data, and the corresponding simulated values, Hy m f : [19] K denotes the Kalman gain and H the measurement operator. D = (d 1 , .., d N ) represents the data matrix holding the data vectors, d m , constructed from the assimilated data d and the error of the assimilated data by assuming that d m are normal distributed about d, according to the error covariance of the assimilated data, R. Here, d m includes the C anth data from the different oceanographic sections and Hy m f , the corresponding C anth values simulated with the Bern3D model for the prescribed set of 22 scaling parameters (y m f ) and for ensemble member (m). The error of the C anth data are assumed to be independent. The vectors d m of the matrix D are computed from d using a random number generator assuming a data error of ± 10%. This error is in-between the error assumed by Sabine et al. [2004] and Matsumoto and Gruber [2005].
[20] The Kalman gain K is computed from the error covariance of the ensemble and the error covariance of the assimilated data R. An estimate of the ensemble covariance, P f , is obtained by where y f is the ensemble mean state, here taken to be as an estimate of the ''true'' model state and superscript T indicates that the matrix is transposed.
[21] A solution for the Kalman gain can be written following Evensen [2003]: [22] In practice, the Kalman gain is computed using the Bern3D model output data using the algorithm described by Evensen [2003], avoiding an explicit quantification of the measurement operator H. In each simulation, air-sea fluxes of C anth are prescribed for 22 regions over the whole ocean. The regions are shown in Figure 1 and Table 1. The prescribed C anth flux is computed for each region as the product of a spatiotemporal pattern describing seasonal variability, a scaling function that accounts for the time evolution of atmospheric CO 2 , and the scaling parameter that is optimized with EnKF and that determines the magnitude of the flux.
[23] The seasonal pattern, P l (i, j, t s ), is based on a CO 2 airsea flux climatology of Takahashi et al. [2002]. The indices, i and j denote the x and y location and time, t s denotes the time of the year. The CO 2 air-sea flux climatology of Takahashi is normalized for each region to yield a unit flux over 1 year. The temporal evolution of the flux of C anth , f(t) is assumed to be proportional to the atmospheric CO 2 increase [Gloor et al., 2003]: t 0 and t 1 denote the start (1770 AD) and the end (2005 AD)of each Bern3D simulation.
[24] The injected flux for region l for the ensemble member m, F m,l , then reads F m;l ði; j; tÞ ¼ P l ði; j; t s ÞfðtÞy m ðlÞ ð 5Þ [25] The initial values of the whole ensemble of scaling parameters for region l, y(l) have been assigned to be normally distributed around zero. A test simulation with 64 instead 32 ensemble members yield deviations from the standard air-sea fluxes of less than 0.02 GtC a À1 except in the two southern subpolar regions, where deviations are up to 0.09 GtC a À1 .

Bern3D Ocean Model
[26] The Bern3D ocean model [Müller et al., 2006] is a cost-efficient three dimensional global circulation model based on the ocean model of Edwards and Marsh [2005]. It has a horizontal resolution of 36 Â 36 grid boxes and a vertical resolution of 32 layers, which are logarithmically spaced. The upper most layer has a thickness of 39 meters, the deepest of around 397 meters. The model is tuned toward data-based CFC and radiocarbon inventories. Model results are found to be in good agreement with observed distribution of different tracers [Müller et al., 2006Parekh et al., 2008;Tschumi et al., 2008]. The model is forced with seasonal fields for temperature and salinity, as well as for the wind stress. The model is able to represent the broad-scale ocean circulation, despite its relatively coarse resolution. Shortcomings of the standard model setup include a too shallow penetration of North Atlantic Deep Water (NADW), a weak gyre circulation and a weak Antarctic Circumpolar Current (ACC). The model's adequate representation of the physical transport timescales of vertical exchange and the model's fast integration time makes it well suited for ensemble simulations.

Bern3D Model Configurations
[27] We apply the Bern3D in different model configurations to investigate uncertainties connected with model transport. The different model setups have been tailored to address known shortcomings in the circulation of the Bern3D ocean model. First, the standard setup, is the model version described by Müller et al. [2006]. Second, the ACC Â 3 setup is the standard version modified by increasing the ACC by a factor of three and by introducing a salt flux from the Pacific to the Atlantic. Increasing the strengths of the ACC leads also to a northward expansion of the ACC and a weakening of the southern subtropical gyres in all basins. The introduction of a salt flux corresponds to a freshwater flux from the Atlantic to the Pacific and increases the formation and propagation of NADW. Third, the PSI Â 3 setup is the standard setup modified by increasing the barotropic streamfunction globally by a factor of three and applying the same salt flux from the Pacific to the Atlantic as in the ACC Â 3 setup. This configuration is characterized by strong vertical and horizontal mixing of tracers and a North Atlantic Deep Water formation of 22 Sv. Fourth, the High Diffusion setup is the standard setup, but diapycnal diffusivity is increased by a factor four from 1 Â 10 À5 ms À2 to 4 Â 10 À5 ms À2 for all passive tracers; the advection and convection for the tracers is the same as in the standard setup. The standard setup yields the best agreement between the modeled and observed distribution of CFC-11 (Table 2).

Results
3.1. Regional Air-Sea Fluxes of C anth : Results for the Global Ocean and the DC*-Based Reconstruction [28] In this subsection, the pattern of air-sea fluxes of anthropogenic carbon as inferred from the DC* reconstruction from the GLODAP database is discussed and compared with the results of Mikaloff Fletcher et al. [2006], who used the same input data in their inversion (Figure 2). The EnKF yields a global uptake flux of anthropogenic carbon of 1.95 GtC a À1 for the year 1995. A strong uptake of anthropogenic carbon is found for the Southern Ocean, particularly in the subpolar regions. The EnKF yields also considerable uptake rates in the midlatitudes of the Indian Ocean. In the Atlantic, most of the C anth enters the sea in northern and midlatitudes. The air-sea flux in midlatitudes and tropical regions for all basins is estimated to be 1.06 GtC a À1 , the uptake flux into the southern polar and subpolar ocean is 0.72 GtC a À1 , and the flux into northern high-latitude regions is 0.17 GtC a À1 . The results of the EnKF compare well with the results of Mikaloff Fletcher et al. [2006] ( Figure 2). In most regions, the two studies agree within The four model configurations are described in section 2.2.3. The results inferred from the TTD, CFC, and TrOCA C anth data sets are obtained with the Standard setup of the Bern3D model. The results with the four Bern3D circulation setups are inferred from the GLODAP [Key et al., 2004] DC* C anth data. Best estimates and their uncertainties (± 1 standard deviation) are shown under ''this study''. Best estimates are calculated from the TTD, CFC, and TrOCA and DC* results for the Standard EnKF-Bern3D setup and using the skill scores of each inversion as weights. The standard deviations from the four inversions (calculated with the skill scores) reflect uncertainties in the C anth data only and are listed in parenthesis. They are added by Gaussian error propagation to the uncertainties from ocean model transport of Mikaloff Fletcher et al. [2006]. In comparison the values from the ocean inversion study of Mikaloff Fletcher et al. [2006] with its standard deviation are also shown. The RMSE of the oceanic distribution of anthropogenic carbon and CFC-11 are also shown. The values are in mmol l À1 and pmol l À1 , respectively. The RMSE of CFC-11 was calculated with data from GLODAP [Key et al., 2004] from the years 1991 -1996. The RMSE C anth was derived from all data points used to constrain the optimization. the uncertainties from ocean transport and agreement is especially good in the Atlantic Ocean.
[29] Four different model setups of the Bern3D model permit us to probe the sensitivity of results to ocean transport ( Figure 2 and Table 2). Although, the root mean square error (RMSE) for CFC-11 and C anth are similar for the four setups, the inferred regional fluxes of C anth are not necessarily in good agreement. Mikaloff Fletcher et al. [2006] has quantified uncertainties related to ocean transport by using transport functions from a suite of ocean circulation models that feature quite different circulation, mixing, and convective regimes. Interestingly, the ranges in inferred regional fluxes for the four Bern3D versions are very similar to the uncertainty range provided by Mikaloff Fletcher et al. [2006]. This suggests that the four different Bern3D setups roughly cover the uncertainty range in ocean transport. Next, deviations between the observation-based [Key et al., 2004] and modeled distribution of CFC-11 are compared (Table 2). Typical CFC-11 concentration range from 50 to a few tenth pmol l À1 in the top 1000 meter of the ocean. Generally good agreement between data and observations is found in forward simulations with the Bern3D model [see also Müller et al., 2006, Figures 1, 10, and 14]. The RMSE is within the range of 0.8 to 0.9 pmol l À1 for the setups. The standard setup yields the smallest and the High Diffusion setup the largest RMSE of the four model versions and for CFC-11. The RMSE between the optimized distribution of C anth and the assimilated C anth data is in the range of 7.11 to 7.75 mmol l À1 . This is comparable to the uncertainty of the assimilated C anth data.
[30] The spatial distribution of the deviations between optimized and assimilated C anth data provides further evidence for the quality of the EnKF and points also to a few weaknesses. For the standard setup, zonal-mean deviations ( Figure 3) are small between 60°S and 40°S and between 20°N and 60°N. Relatively large deviations are found in the upper thermocline around the equator and around 25°S between 350 and 600 m. We tentatively attribute this to a too weak formation and extension of Antarctic Intermediate Water. In the setup with increased barotropic streamfunction, horizontal exchange is enhanced, the equatorial deviations in optimized and assimilated C anth reduced, and the air-sea flux in the subpolar and polar regions increased by 0.24 GtC a À1 relative to the standard setup ( Table 2). The deviations between model results and assimilated data are also reduced in the southern polar region, although the airsea flux is higher than for the standard case. Taken together, the inferred uptake in southern subpolar and polar regions might be somewhat biased on the low bound for the standard model version. The air-sea flux in midlatitude regions is increased and the flux in the Southern Ocean reduced relative to the standard in the version with a high diapycnal diffusivity for passive tracers. This also yields a reduction in the deviations between modeled and assimilated data in the tropical thermocline. However, this solution is not considered to be particularly realistic, as the thermocline appears to be too diffusive. Turning to the North Atlantic, the four different model setups yield only relatively small difference in the EnKF optimization. This suggests that the applied modification in the freshwater flux between Atlantic and Pacific have little impact on inferred air-sea fluxes of C anth , although the extension of NADW and strengths of NADW formation is more realistic with an enhanced freshwater export from the North Atlantic.
[31] In summary, the combination of the EnKF and the Bern3D model is able to reproduce the observation-based C anth distribution as reconstructed with the DC* method reasonably well. The four Bern3D model versions provide a quantification of uncertainties related to ocean transport in C anth assimilations and uncertainty ranges that are similar to those obtained with a suite of circulation models. The global air-sea flux of C anth in 1995 is estimated to be 1.95 to 2.16 GtC a À1 across the model range. The optimized uptake flux in polar and subpolar regions of the Southern Ocean is in the range of 0.45 to 0.96 GtC a À1 for the four setups and the higher estimates are considered to be more realistic. Finally, inferred regional air-sea fluxes for the Atlantic are robust across model versions and agree favorably with the results from Mikaloff Fletcher et al. [2006].

Regional Air-Sea Fluxes, Storage, and Meridional Transport of Anthropogenic Carbon 3.2.1. Results for the Global Ocean and Four Different C anth Reconstructions
[32] Uncertainties arising from uncertainties in the reconstruction of C anth are addressed here. The C anth fields from the TTD, CFC, and TrOCA methods are assimilated in the Bern3D standard setup and the results are compared to those from the DC* method (Figure 2 and Table 2). There is broad agreement between the results from the four inversions. However, inferred air-sea fluxes, changes in storage, and transport divergence of C anth deviate notably between methods in various regions. The deviation between assimilated data and optimized results provides a check on the quality of the assimilation (Table 2 and Figure 3). Deviations, expressed as RMSE, can arise both from shortcomings of the ocean transport model as well as from inconsistencies and biases of the reconstructions. The small RMSE between modeled CFC-11 and observations is encouraging and suggests a relatively high quality of the model transport field. We note that the model has also been tested with a suite of other transient and steady state tracers [Müller et al., 2006Tschumi et al., 2008;Parekh et al., 2008]. Model related shortcomings are, as discussed in the previous section, a weak penetration of North Atlantic Deep Water and weak formation of intermediate waters. The TTD-based reconstruction yields the smallest RMSE (6.8 mmol l À1 ). The RMSE for the DC* reconstruction is 7% larger, and those for the TrOCA and CFC reconstructions are 45% and 100% larger, respectively. The RMSE for the DC* reconstruction corresponds well to the estimated precision of 10 m mol kg À1 of the method [Gruber et al., 1996]. As indicated earlier, the basic assumption of a purely advective flow (discrete transit time) underlying the CFC reconstruction does not appear realistic. The CFC reconstruction yields substantial C anth concentrations in the deep Indian and Pacific, resulting in zonal-mean residuals between data and model of around 15 mmol l À1 at depth.The RMSE for the TrOCA method is 70% larger than the estimated uncertainty of the method.
[33] Skill scores are used as weights for the individual inversion results to compute best estimates and their uncertainties (Tables 2 and 3). The variance, together with the correlation between optimized and reconstructed fields, are used to compute skill scores, S, for each inversion following Taylor [2001]: where s is the normalized variance and R the correlation between modeled and data-based distribution of C anth . The maximum correlation attainable, R 0 , is set to 1. Global C anth uptake for 1995 is in the range of 1.95 to 2.30 GtC a À1 , with a best estimate of 2.10 ± 0.16 GtC a À1 . The best estimate is 0.85 ± 0.17 GtC a À1 (range: 0.72 to 1.09 GtC a À1 ) for the polar and subpolar Southern Ocean, 1.07 ± 0.18 GtC a À1 (0.84 to 1.17 GtC a À1 ) for the tropical and midlatitude regions, and 0.18 ± 0.06 GtC a À1 (0.12 to 0.27 GtC a À1 ) for the northern high latitude. Interestingly, deviations between TTD and DC*-based results are equal or smaller than 0.11 GtC a À1 for these aggregated areas. This is surprising, as column inventory for the DC*-based method are much smaller than for the TTD method in the Southern Ocean ( Figure 1).
[34] Turning to individual regions, we note that the largest spread in results is found for the Southern Ocean regions (Figure 2a). The partitioning in uptake among the subpolar and polar regions does not appear well constrained. For example, uptake in the relatively small subpolar Atlantic region is 0.12 GtC a À1 smaller for the TTD reconstruction than for the DC* reconstruction, whereas uptake in the polar Southern Ocean is higher by about the same amount for TTD than for DC*. This suggests that current estimates of uptake fluxes in individual Southern Ocean regions are uncertain by up to a factor of two. A different allocation of flux is also evident in the Pacific. The TTD reconstruction yields an uptake in the tropical Pacific that is 70% (0.15 GtC a À1 ) higher than that from the DC* reconstruction, whereas inferred uptake in the southern midlatitude Pacific is almost a factor of three lower for TTD than for DC*.
[35] The optimized changes in C anth storage (Figure 2b) agree well and differences among methods are smaller than 0.06 GtC a À1 for individual regions for the nominal year 1995. Turning to individual basins, storage rates are within 0.60 (TTD) to 0.70 (TrOCA) GtC a À1 for the Atlantic and 0.30 (TrOCA) to 0.36 (CFC) GtC a À1 for the Indian. For the Pacific, the inferred change in storage is lower for the DC* reconstruction (0.65 GtC a À1 ) than for the TrOCA and TTD (0.73 and 0.74 GtC a À1 ) and the CFC (0.83 GtC a À1 ) reconstructions. A characteristic feature of all inversions is the large export of C anth from the subpolar regions into southern midlatitude regions and toward the tropics (Figure 2c). Only about 40% (range 34 to 54%) of the uptake from the atmosphere remains in the subpolar region. Meridional transport of anthropogenic carbon (Table 4) in the Atlantic is northward for the nominal year 1995 and broadly consistent across all reconstructions. Meridional transport is converging for most Atlantic regions, except  The values in parentheses in the DC* column indicates values obtained by constraining the inversion with all available data from the GLODAP data set, and not only with data from the four sections. The average (all methods) is the mean of all six methods with the standard deviation. The average (four methods) column is the average from the DC*, 8 À C T 0, TDD, and TrOCA reconstructions.
in the subpolar Atlantic and in the southern tropical Atlantic, two regions with high uptake from the atmosphere. Meridional transport convergence is high in the northern low-latitude Atlantic. The northern low-latitude, midlatitude, and high-latitude regions show the largest change in carbon storage of all Atlantic regions as a result of high transport convergence and uptake from the atmosphere. In the Pacific, we find that C anth is transported from the equator toward the midgyre regions and that transport convergence and storage rate are high in southern midlatitudes.
[36] Differences in inferred regional transport divergences ( Figure 2c) and meridional transport rates (Table 4) are substantial in some regions across methods. Largest differences are again found in the Southern Ocean. The polar Southern Ocean is losing some C anth by northward transport in the DC*-based inversion, whereas the TTD-based inversion suggest an influx of 0.1 GtC a À1 . Similar differences between the two methods are found for the subpolar Atlantic and the southern midlatitude Pacific. The TrOCA derived divergences and air-sea fluxes appear anomalous for the subpolar Indian and Pacific and the southern midlatitude Indian and Atlantic when compared to results from the other three reconstructions. For example, air-sea fluxes and divergences are more than 0.2 GtC a À1 higher for TrOCA than any of the other three methods in the subpolar Indian and Pacific region. The patterns of air-sea flux, meridional transport, and storage of C anth , as inferred with the EnKF, are consistent with oceanographic understanding. For example, the large C anth uptake and subsequent northward transport in the southern subpolar region is linked to the northward spreading of Antarctic Intermediate Water. Relatively large uptake of anthropogenic carbon are also obtained for the northern midlatitude and high-latitude Atlantic regions (46°N to 71°N) and in the southern tropical Atlantic (3°S to 19°S). These fluxes are related to the formation of NADW and convective activity in high latitudes, as well as the upwelling of anthropogenic carbon depleted waters in the tropics. Uptake rates in the midgyre regions of the north and south Atlantic are low.
[37] Meridional transport from all reconstructions is broadly consistent with the rates estimated by Mikaloff Fletcher et al. [2006] and hydrographic data-based estimates [Rosón et al., 2003;Macdonald et al., 2003]. In agreement with Mikaloff Fletcher et al. [2006], we find a much higher C anth uptake by air-sea exchange north of 25°N in the Atlantic than Macdonald et al. [2003], where the air-sea flux is estimated to be only 0.02 GtC a À1 .
[38] In summary, results for the TTD and DC*-based inversion, the two inversions yielding the lowest RMSE, agree well on large aggregated regions. However, differences are substantial in individual regions such as in the Southern Ocean and the tropical Pacific.

Results for the Atlantic and for Six Different C anth Reconstructions
[39] C anth reconstructions from the IPSL [Lo Monaco et al., 2005] and the 8 À C T 0 are available on a few selected Atlantic transects in addition to reconstructions with the TTD, CFC, TrOCA, and DC* method. These section data have been applied in the EnKF. We will first discuss reliability of these ''Atlantic-only'' inversions, before turning to individual results and their spread across reconstruction methods.
[40] Boundary conditions have to be prescribed for the Atlantic-only inversion. The IPSL and 8 À C T 0 results are only available for selected sections (Figure 1). The air-sea fluxes for the polar Southern Ocean and the different regions of the Indian and Pacific were prescribed using the results from the DC*-based inversion; these DC*-based air-sea fluxes are prescribed for all six methods for consistency. Limited data coverage and boundary conditions could affect results. The comparison between results from the global DC* inversion discussed previously (Table 2) and the Atlantic-only DC* inversion (Table 3) reveals that deviations in regional air-sea fluxes of C anth are relatively small ( 0.04 GtC a À1 ). Deviations are also small for the CFC, TTD, and TrOCA global and Atlantic-only inversions, the southern midlatitude and the subpolar Atlantic regions excepted. In these two regions, results for the CFC, TTD, and TrOCA reconstruction appear affected by the DC*-derived boundary conditions prescribed in the Atlantic-only inversion. Uptake fluxes are lower in southern midlatitudes and higher in the subpolar Atlantic in the Atlantic-only inversion compared to the global inversions. Increasing the prescribed air-sea flux for the polar Southern Ocean by a factor of three reduces uptake in the subpolar Atlantic, but air-sea fluxes and their standard deviation from the range of reconstructions remains unaffected in other regions.
[41] We conclude that the data from the four Atlantic sections permit us to analyze the spread in results arising from differences in the six C anth reconstructions, though the magnitude of individual results is somewhat affected by limited data coverage and prescribed boundary fluxes. The RMSE between the assimilated and optimized Atlantic C anth data (Table 3) is (almost 13 mmol l À1 ) quite large for the IPSL and CFC-shortcut reconstructions, whereas the range in the RMSE is 7.2 to 8.7 mmol l À1 for the other four reconstruction methods. The best agreement between assimilated and optimized model distribution is obtained for the TTD-based C anth data.
[42] The high RMSE for the CFC-shortcut and the IPSL reconstructions arise for different reasons. The CFCshortcut C anth data decline on average much more smoothly toward the deep ocean than those from the other reconstructions. The model is not able to reproduce the large penetration depth of C anth suggested by the CFC-shortcut method as this would require a vigorous surface-to-deep mixing. The IPSL C anth data suggest high concentrations at the surface and a strong decline with increasing depth. This spatial pattern is again not reproducible by the model, and results in a high RMSE. In the following, we will focus discussion on the reconstructions with a low RMSE. There is a large spread in the total and regional air-sea fluxes of C anth in the Atlantic inferred from the six different C anth reconstructions ( Figure 4 and Table 3). The optimized Atlantic uptake in 1995 is in the range from 0.54 to 0.79 GtC a À1 for all six methods and from 0.54 to 0.66 GtC a À1 , when omitting the CFC and IPSL results. The TTD and the 8 À C T 0 C anth data provide the lower bound, while the IPSL C anth reconstruction yields the upper bound and the TrOCA, DC*, and CFC-shortcut methods yield intermediate estimates of the total Atlantic uptake.
[43] The DC* and 8 À C T 0 C anth data EnKF assimilation suggests that 29% of the total Atlantic uptake is in the subpolar box (63°S to 43°S). In contrast, the other four C anth reconstructions yield a contribution of about 43% from the subpolar region. Not only the relative contribution, but also the absolute flux for the subpolar region shows a large spread across methods. Large relative deviations in uptake of up to a factor of five (factor three without IPSL and CFC) are found for the southern tropical region. Correspondingly, the spread in transport divergence is from 0.01 to 0.18 (0.01 to 0.11 without IPSL) GtC a À1 and 20 to 70% of the carbon taken up from the atmosphere is exported by meridional transport. Interestingly, the ranking in total uptake in the Atlantic for the different reconstructions is not as expected on the basis of the magnitudes of the column inventories. For example, the CFC-shortcut method yields by far the largest column inventories (Figure 1), but the Atlantic uptake is lower than those derived from the IPSL and DC* C anth data. This demonstrates that not only the total inventory of anthropogenic carbon, but also its distribution within the ocean codetermines the inferred magnitudes of the air-sea flux. This suggests that it might be difficult to assess the uncertainties arising from uncertainties in C anth by conventional sensitivity and error propagation analyses. Analyzing results from different C anth reconstructions as done here seems more reliable to estimate uncertainties.
[44] In summary, uncertainties in the determination of the concentration of C anth in the ocean yield a large spread in the EnKF inferred air-sea fluxes and meridional transport rates of anthropogenic carbon. The uncertainty range in optimized Atlantic air-sea flux is 20% of the average Atlantic uptake for the DC*, TTD, TrOCA and 8 À C T 0 reconstructions.

Discussion and Conclusion
[45] A key challenge of any scientific analysis is not only to provide best estimates, but to assess and quantify uncertainties of results. Considering different sources of uncertainty usually results in larger error bars and a less well quantified situation. This is not particularly appealing, but necessary to avoid false interpretations of results. In this work, we contribute to this task of quantifying uncertainties. An Ensemble Kalman Filter has been applied to assimilate anthropogenic carbon data from four different global and six different Atlantic C anth reconstructions into the Bern3D dynamic ocean model. Uncertainties in inferred regional airsea fluxes, meridional transport, and ocean transport divergences of anthropogenic carbon arising from uncertainties in the ocean distribution of C anth have been explored in detail. Uncertainties arising from uncertainties in ocean transport have been discussed. The result is a new set of best estimates for air-sea and meridional fluxes of anthropogenic carbon with more comprehensive estimates of uncertainties than provided by previous assessments. There is broad agreement between the results from inversions of the different C anth reconstructions. However, inferred airsea fluxes and transport fluxes of C anth within the ocean deviate notably between methods in various regions. Global uptake of C anth is estimated to be 2.10 ± 0.30 GtC a À1 for the nominal year 1995 and 131 ± 18 GtC over the period 1770 to 2000 AD. The largest uptake is found in the subpolar regions of the Southern Ocean. Consistent with oceanographic understanding, most of the anthropogenic carbon that entered the subpolar region by air-sea flux is transported northward and carbon taken up in the equatorial Pacific is transported toward midgyre regions.
[46] The largest spread in results is found for the Southern Ocean regions. The partitioning in uptake between the subpolar and polar region is not well constrained. Uptake fluxes in individual Southern Ocean regions are uncertain by up to a factor of two. Uncertainties are also substantial in the Pacific. For example, the TTD reconstruction yields a C anth uptake in the tropical Pacific that is 70% higher than that from the DC* reconstruction, whereas inferred uptake in southern midlatitude Pacific region is almost a factor of three lower for the TTD than for the DC* reconstruction.
The results of any inversion are only as reliable as the input data. All reconstructions assume that ocean transport remained in steady-state over the industrial period. The DC*, IPSL, TrOCA and 8 À C T 0 also assume that preformed concentrations of a range of tracers remained invariant over time (the effects of global warming and of estimated changes in CaCO 3 production is included in the 8 À C T 0 reconstruction (M. Vázquez-Rodríguez et al., manuscript under review, submitted 2008)). The impact of these assumptions remains yet to be quantified.
[47] There are methodological uncertainties and potential biases that are specific to individual reconstructions. There is a rich literature available discussing the benefits and shortcomings of individual methods. Here, we just briefly highlight a few findings without intention for completeness. Different explicit or implicit assumptions on ocean transport and mixing and how the preformed tracer concentrations and air-sea CO 2 disequilibria are computed in detail are a source of difference. For example, the DC* method assumes that transport and mixing occurs along isopycnals [Gruber et al., 1996]. Six end-members were defined to resolve mixing in the East Atlantic in the 8 À C T 0 reconstruction (M. Vázquez-Rodríguez et al., manuscript under review, submitted 2008) and a single dominant water source is assumed in the TTD and CFC methods [Waugh et al., 2006]. Comparison of C anth data from the 8 À C T 0 method and the DC* reveal that such differences cause large discrepancies in the Southern Ocean C anth inventories (see, e.g., Figure 1). Three of the methods (TTD, DC*, and CFC-shortcut) have been assessed within the framework of an ocean general circulation model [Matsumoto and Gruber, 2005;Waugh et al., 2006;Matear et al., 2003]. Waugh et al. [2006] find that estimates of C anth from the TTD method are biased toward the higher end in the Southern Ocean because Figure 4. (a) Atlantic air-sea fluxes, (b) storage rates, and (c) transport divergence (uptake minus storage) in GtC a À1 scaled to the year 1995 derived from six different data-based estimates. The star on the DC* column in Figure 4a shows results obtained by assimilating all available DC*-based C anth data in the Atlantic. The error bar in Figure s 4a, 4b, and 4c is one standard deviation around the average from the different methods.
of the assumption of a temporally constant air-sea CO 2 disequilibrium and that predicted C anth closely matches the directly simulated distribution in other areas. Matsumoto and Gruber [2005] find that the DC* method tends to overestimate C anth in relatively young water, but underestimates it in old water. In the global mean, a positive bias of 7% is estimated, but local biases can be larger. Main sources of these biases are the assumption of a time invariant air-sea disequilibrium, errors in identifying different end-member water types, and the application of a single water mass age to determine the air-sea disequilibrium in young waters. The latter leads to a discontinuity in vertical profiles of C anth [Waugh et al., 2006]. The analysis of Matear et al. [2003] and Waugh et al. [2006] suggests that the CFC-shortcut method is not reliable for old waters such as those found in the deep and in the Southern Ocean. The IPSL, TrOCA and 8 À C T 0 reconstruction have not yet been evaluated in the same stringent manner.
[48] The different EnKF inversions of C anth are not of equal quality. This is indicated by different root mean square errors and skill scores computed from the deviations between assimilated data and optimized model fields (Table 2). A large RMSE may either be caused by deficiencies in the assimilation scheme, the ocean transport model or by problems in the underlying C anth reconstruction. The RMSE (7 mmol l À1 ) are reasonably small for the global inversion of the TTD and DC* C anth reconstructions, but almost twice as large as the inversion with the CFC data. The regional air-sea fluxes inferred with the EnKF-Bern3D approach from the GLODAP DC* data agree within uncertainties with the estimates of Mikaloff Fletcher et al. [2006] who assimilated the same C anth data into a range of models using a Green's Function approach. The Bern3D model also reproduces the distribution of CFC-11 well. We conclude that both the Green's Function approach and the Ensemble Kalman Filter are well suited to assimilate observationbased data into ocean models and that the quality of the Bern3D tracer transport is comparable to those of the general ocean circulation models used by Mikaloff Fletcher et al. [2006].
[49] Interestingly, the ranking in total uptake in the Atlantic for the different reconstructions is not as expected from the magnitudes of the column inventories. For example, the TTD reconstruction yields a similar or even higher C anth column inventory than the DC* C anth reconstruction, but has a lower Atlantic uptake than the latter one. This demonstrates that not only the inventories, but also the gradients within the ocean determine the inferred fluxes. An implication is that conventional Gaussian error propagation to estimate the uncertainty in C anth data from uncertainties in Redfield ratios and concentration data may not recover the full uncertainty range associated with uncertainties in C anth data. Arguably, a better approach to estimate uncertainties from uncertainties in input data is to use data from different reconstruction methods, each with its own set of assumptions and potential systematic biases. Mikaloff Fletcher et al. [2006] concluded that inverse flux estimates generally tend to be more sensitive to the choice of model than to biases in the anthropogenic CO 2 estimates. However, assimilating data from six different published reconstruction methods of C anth yields a spread in inferred air-sea fluxes comparable to or larger in magnitude than that from different models and model setups. The standard deviations in inferred Atlantic air-sea fluxes from different models and model setups given by Mikaloff Fletcher et al. [2006] are in the range of ±0.01 to ±0.04 GtC a À1 for the Atlantic regions (regions 1 to 7), whereas the standard deviations from the six different reconstruction methods vary from ±0.01 to ±0.07 GtC a À1 (±0.01 to ±0.05 GtC a À1 without IPSL and CFC) ( Table 3). We conclude that the uncertainty estimates by Mikaloff Fletcher et al. [2006] might be biased somewhat in the low end and that both uncertainties in reconstructed C anth and uncertainties in ocean model transport must be included in any error assessment.
[50] In conclusion, the uncertainties in the C anth reconstructions indicates that we cannot yet be confident of how much anthropogenic carbon is entering, stored in, and exported from the Southern Ocean. The Southern Ocean is important for projecting future atmospheric CO 2 and global warming, yet forward ocean model results and inverse analyses show the largest discrepancies in Southern Ocean regions.