Interactive comment on “ A Joint Global Carbon Inversion System Using Both CO 2 and CO 2 Atmospheric Concentration Data

This manuscript presents an interesting approach for using 13CO2 data as extra constraints for top-down flux inversions based on in-situ surface CO2 data. This approach has taken into account spatial variation of isotropic discrimination and disequilibrium by using a terrestrial biosphere model and an ocean model to simulate the discrimination rates. The manuscript is well written, and their results are interesting. It should be published after minor revision.


Introduction
Over the last few decades, much progress has been made in estimating the global carbon cycle using different methods (Houghton et al., 2007;Canadell et al., 2007;Le Quéré et al., 2013).In particular, atmospheric CO 2 mole fractions measured near the surface have been used to infer the carbon flux over land and ocean surfaces through atmospheric inversion (Rödenbeck et al., 2003;Michalak et al., 2005;Peylin et al., 2005;Peters et al., 2007).However, the uncertainty in the inferred flux is still very large, mostly because of the insufficient number of observation stations and the error in modelling the atmospheric transport of CO 2 from the surface to the observation stations.To reduce this uncertainty, it would be useful to introduce constraints to the inversion using other gas species that are associated the CO 2 flux.
Measurements of the atmospheric concentration of the stable isotope 13 CO 2 at a number of stations across the globe since 1994 have been compiled in a database (GLOBALVIEW-CO2C13, 2009), and the number of extended 13 CO 2 records from January 1994 to January 2009 increased to 76 by 2009.The mole fraction of 13 CO 2 to CO 2 in the atmosphere is about 1.1 %, and the CO 2 exchange between the surface and the atmosphere generally induces concurrent 13 CO 2 exchange.However, the proportion of the 13 CO 2 flux relative to the CO 2 flux differs at dif-Published by Copernicus Publications on behalf of the European Geosciences Union.
ferent locations and different times due to different mechanisms that discriminate against heavier 13 CO 2 molecules in the exchange processes, and therefore the 13 CO 2 concentration measured in the atmosphere contains additional information for the CO 2 flux.This information is useful for differentiating between terrestrial and oceanic CO 2 exchanges with the atmosphere because the terrestrial CO 2 flux experiences much greater discrimination against 13 CO 2 than does the oceanic CO 2 flux (Tans et al., 1990;Ciais et al., 1995a;Francey et al., 1995).Observed 13 CO 2 mole fractions can also provide independent information on the net CO 2 exchange over land and ocean because the net carbon flux to the surface discriminates against heavier 13 CO 2 (Fung et al., 1997;Randerson et al., 2002;Suits et al., 2005).The 13 CO 2 observations over the globe, albeit with a limited number of stations, could therefore be used to assist in quantifying the global carbon cycle.
In previous studies (Siegenthaler and Oeschger, 1987;Keeling et al., 1989a;Francey et al., 1995;Randerson et al., 2002), atmospheric 13 CO 2 observations have been used to separate ocean and land CO 2 fluxes through the use of a technique dubbed "double deconvolution", by which the CO 2 fluxes of land and ocean are separated (deconvolved) based on different discrimination rates against 13 CO 2 in the atmospheric CO 2 exchange with land and ocean surfaces.This double deconvolution often assumes that the discrimination rates over land and ocean are spatially uniform, although they can be temporally variable.Through forward atmospheric transport modelling, the ocean and land CO 2 fluxes were also separated based on the spatial gradients of the measured 13 CO 2 / CO 2 ratio either globally (Keeling et al., 1989b) or by latitudinal bands (Ciais et al., 1995a).The same 13 CO 2 data have also been used in inverse modelling of the surface CO 2 flux (Enting et al., 1995;Rayner et al., 1999Rayner et al., , 2008)).Enting et al. (1995) pioneered a methodology for inverting annual mean ocean and land CO 2 fluxes from both atmospheric CO 2 and 13 CO 2 concentration data for 12 ocean regions and 8 land ecosystems for the 1986-1987and 1989-1990periods. Rayner et al. (1999) ) developed a different methodology to invert monthly CO 2 fluxes for 12 ocean and 14 land regions for the period from 1980 to 1995 from CO 2 observations at 12 stations and 13 CO 2 and O 2 / N 2 observations at 1 station.Rayner et al. (2008) refined their methodology and applied it to the period from 1992 to 2005 using CO 2 at 67 sites and 13 CO 2 at 10 sites.These studies showed the usefulness of the additional information from 13 CO 2 observations in improving the inversion of annual mean and seasonality of the CO 2 flux over land and ocean.In these inversion studies, the discrimination rate for land is either assumed to be a constant (Enting et al., 1995;Rayner et al., 1999) or allowed to vary with the areal fraction of C4 plant in a region (Rayner et al., 2008).These inversions based on the Bayesian principle were also constrained with only simple prior estimates of the terrestrial and oceanic CO 2 and 13 CO 2 fluxes.Since the data density (the numbers of CO 2 and 13 CO 2 observation sites) is low, the assumed discrimination constants and these prior estimates would have considerable influence on the inverted results, as this is clearly demonstrated in Enting et al. (1995).
Atmospheric CO 2 observations have been extensively used to estimate the carbon flux over ocean and land through inverse modelling using Bayesian synthesis (Gurney et al., 2002;Rödenbeck et al., 2003;Baker et al., 2006;Peylin et al., 2005) or data assimilation techniques (Peters et al., 2007;Zhang et al., 2014).Atmospheric inversion studies (Gurney et al., 2003;Jacobson et al., 2007) often produced ocean sinks considerably smaller than those estimated based on observed gradients in dissolved inorganic carbon (DIC) in interior ocean using ocean circulation models (Steinkamp and Gruber, 2013).Recent estimates for the ocean sink for anthropogenic CO 2 in the 2000s were based on DIC ranges from 1.6 to 2.6 Pg C year −1 (Park et al., 2010;Wanninkhof et al., 2013;Landschützer et al., 2014;Majkut et al., 2014;DeVries, 2014;Rödenbeck et al., 2014) with an uncertainty of about 0.6 Pg C year −1 , while atmospheric inversion results are not yet reliable enough to be included in a global ocean sink synthesis (Le Quéré et al., 2013).The partition between ocean and land fluxes using atmospheric inversion techniques is sensitive to errors in atmospheric transport modelling (Patra et al., 2005;Baker et al., 2006;Stephens et al., 2007) and prior fluxes for land and ocean used to constrain the inversion (Zhang et al., 2014;Chen et al., 2015).It would therefore be highly desirable to use 13 CO 2 observations to constrain this partition in the inversion process.An accurate partition between ocean and land sinks is important in global carbon cycle research because (1) land sinks are still more reliably estimated as the residual of the global carbon budget than those from land-based data (Le Quéré et al., 2013) and (2) ocean sink estimates based on DIC in ocean water also suffer from considerable errors due to insufficient DIC observations and ocean circulation modelling (DeVries, 2014).
The overall goal of this study is to explore the information content of 13 CO 2 measurements for global CO 2 flux estimation through developing a Bayesian synthesis inversion system that uses both CO 2 and 13 CO 2 observations.This system is effectively a new double deconvolution system with the capacity of considering the spatial variations of the prior carbon flux and all major isotopic parameters including photosynthetic discrimination, respiratory signature and disequilibrium rate.In this study, this new system is used to achieve the following objectives: (1) to partition between ocean and land sinks with consideration of the spatial distributions of 13 CO 2 isotopic parameters over ocean and land, (2) to evaluate the importance of considering the spatial distributions of the 13 CO 2 discrimination rate over land in the inversion of the CO 2 flux and (3) to assess the impacts of the errors in disequilibrium flux estimation on the flux partition between ocean and land.To achieve these objectives, a terrestrial ecosystem model called the Boreal Ecosystem Productivity Simulator (BEPS) is further developed to simulate the Figure 1.A global nested inversion system with a focus in North America, in which oceans are divided into 11 regions and land areas are divided into 9 large and 30 small regions outside and within North America, respectively.Also shown are CO 2 and 13 CO 2 observation stations included in the GLOBALVIEW database and used in this study; 10 of the stations are marked with their names because they are selected to compare prior and posterior concentrations in Fig. 11.spatial distributions of the 13 CO 2 discrimination and disequilibrium rates over land for use in a global Bayesian synthesis inversion with 13 CO 2 constraint.BEPS is also used to produce CO 2 prior fluxes globally to regularize the inversion.

Inversion system
The nested inversion system with a focus on North America developed by Deng et al. (2007) is adopted in this study.In this system, two of the Transcom regions (Gurney et al., 2002) in North America are divided into 30 regions according to ecosystem types and administrative boundaries (Fig. 1), in order to reduce spatial aggregation errors in the inversion over North America and to investigate the inverted spatial distribution of the carbon flux against ecosystem model results.This nested region serves the purpose of evaluating the influence of the spatial distribution of isotopic discrimination on the inverted carbon flux at a relatively high resolution.Also shown in Fig. 1 are the spatial distributions of 210 CO 2 and 73 13 CO 2 observation sites selected in this study from the NOAA GLOBALVIEW database.Most 13 CO 2 sites (except 11) are collocated with CO 2 sites.

Synthesis Bayesian inversion with CO 2 observations
To estimate the CO 2 flux (f ), we represent the relationship between CO 2 measurements and the flux from the surface by a linear model: where c m×1 is a given vector of m CO 2 concentration observations over space and time (m equals number of stations times number of months, and for CO 2 -only inversion; it is 12 600, i.e. 210 stations × 60 months, [2000][2001][2002][2003][2004]; ε m×1 is a random error vector with a zero mean and a covariance matrix cov(ε) = R m×m ; G m×(n−1) is a matrix representing a transport (observation) operator, where n − 1 is the number of fluxes to be determined (equals 3000, i.e. 50 regions × 60 months, 2000-2004); A m×1 is a unity vector (filled with 1) representing the assumed initial well-mixed atmospheric CO 2 concentrations (c 0 ) before the first month; and f (n−1)×1 is an unknown vector of monthly carbon fluxes of the 50 regions.
Combining matrixes G and A as M m×n = (G, A) and vectors f and c 0 as s n×1 = (f T , c 0 ) T , Eq. ( 1) can be expressed as (2) The inverse problem of estimating s from c is often poorly constrained and a Bayesian approach is used to circumvent this problem.Pre-existing knowledge and models incorporating additional sources of information can be used to provide an initial estimate of s, known as the a priori, to constrain the inversion.This a priori is then updated when it is combined with information from c measurement to form a posterior estimate of s, known as the a posteriori.In Bayesian synthesis inversion (Tarantola, 1987), the following objective function is employed in place of the traditional least-square objective www.geosci-model-dev.net/10/1131/2017/Geosci.Model Dev., 10, 1131-1156, 2017 function: where s p n×1 is the a priori estimate of s, the covariance matrix Q n×n represents the uncertainty in the a priori estimate and R m×m is the transport model-data mismatch error covariance.By minimizing this objective function expressed in Eq. ( 3), we obtain the posterior best estimate of s as (Enting, 2002): Meanwhile the posterior uncertainty matrix for the posterior flux can be deduced as follows: Following the methodology of Deng and Chen (2011), the CO 2 concentration matrix c in the above equations is the residual concentration after subtracting the observed concentration with contributions from fossil-fuel emission, biomass burning, the prior ocean flux and the prior biospheric flux (see Sect. 2.4 for detail).In this way, the values in s p are set to zero and the inverted flux s is considered to be an adjustment to the prior flux that contributes to the pre-subtracted portions of the CO 2 concentration.

Synthesis Bayesian inversion with both CO 2 and 13 CO 2 observations
We attempt to use 13 CO 2 observations to provide an additional constraint to the otherwise CO 2 -only inversion presented above.This additional constraint is possible on the grounds that air 13 CO 2 concentration is affected differently by carbon fluxes from ocean and land surfaces.Since the 13 CO 2 gas is transported passively in similar ways as CO 2 , the same transport matrix M applies to 13 CO 2 data to associate 13 CO 2 observations with the surface 13 CO 2 flux.This simple treatment of the transport matrix differs from Rayner et al. (2008), who considered the reduced response of observed 13 CO 2 concentrations to surface fluxes with time due to its accumulated exchange with the surface.As we are interested in the net CO 2 flux, the exchanges of both 13 CO 2 and CO 2 with the surface are consistently not included in the M matrix calculation, although this simplification would induce errors in the inverted CO 2 flux when the accumulated exchanges are spatially highly heterogeneous.In order to conduct an inversion using both CO 2 and 13 CO 2 observations, we simply append 13 CO 2 -related data to the c, R and M matrixes in Eq. ( 4), while the s matrix remains unchanged as the purpose of this joint inversion is only to optimize the CO 2 flux.For c and R, 13 CO 2 observations and their variances are appended directly to the original matrixes for the CO 2 -only case, as shown in Eq. ( 6).Similarly, the M matrix is also extended to consider 13 CO 2 transport, and the relevant elements for the 13 CO 2 observation stations are from the original M matrix.However, these elements are multiplied by the 13 CO 2 discrimination rate over land or ocean for each region and each month in order to relate the CO 2 flux to the temporal variations in the measured air 13 CO 2 composition at each station and each month.The extended M is a combination of the corrected M matrix appended to the M matrix for CO 2 (see below) where c i is the CO 2 concentration (i = 1 to m) and 13 C composition (i = m+1 to m+k) in the air from the starting month (i = 0); M ij is the transport operator between region-month j (hereafter simply referred as "region") and station-month i (hereafter simply referred as "station"); and W ij = D j M ij , in which D j is the discrimination rate against 13 CO 2 in the CO 2 flux for region j .In the inversion procedure, the difference in concentration between two consecutive times is equated with the flux during the time interval (1 month).
In order to calculate D j and C i (i = m + 1 to m + k) in Eq. 6, some theoretical development is made according to the 13 CO 2 budget equation derived by Tans et al. (1993): , where C a is the CO 2 pool in the atmosphere (in Pg C), δ a is the 13 C composition of the atmosphere in ‰, F f is the carbon emission from fossil fuels and biomass burning, δ f is the 13 C composition of fossil fuels or biomass, F lph is the photosynthetic carbon uptake by the land biosphere (always positive), F lb is the respiratory carbon flux of the land biosphere (always positive), ε lph is the photosynthetic discrimination of the land biosphere in ‰, δ lb is the 13 C composition of the land respiratory carbon flux (see Sect. 2.2.2), δ e lb is the biospheric 13 C composition in equilibrium with the current atmosphere (i.e. in 2003), F oa is the one-way carbon from the ocean surface to the atmosphere (always positive), F ao is the one-way carbon flux from the atmosphere to the ocean surface (always positive), ε ao is the air-ocean fractionation, ε oa is the air-ocean fractionation and δ e a is the 13 C composition in equilibrium with the ocean surface.Equation (7) states that the temporal variation of the measured 13 C composition in the atmospheric CO 2 is determined by contributions from the various sources: fossil fuels and biomass burning (term 1 of the right-hand side of Eq. 7), net land biosphere carbon uptake (term 2), one-way respiratory flux from the land biosphere (term 3), net carbon flux of the ocean (term 4), and one-way ocean-atmosphere flux (term 5).The one-way carbon fluxes from land and ocean surfaces are important sources of 13 C because the atmosphere is in isotopic disequilibrium with these surfaces due to the long-term change of the atmospheric 13 C composition.Similar to other terms in Eq. ( 7), these disequilibrium fluxes are also called isofluxes (Rayner, 2001).
In order to reduce the errors of our inversion system (Eq.6) that assumes linear relationships between fluxes and concentrations, the contributions of all fluxes, including prior biospheric and ocean fluxes, to the CO 2 concentration are subtracted from the measured CO 2 concentration prior to the inversion (Deng and Chen, 2011).Accordingly, the contributions of all 13 C sources to the 13 C concentration in the atmosphere are also subtracted from the measured 13 C concentration.The purpose of the inversion is then to find the residual CO 2 flux, denoted as S in Eq. ( 6).For this purpose, we denote S lN = −(F lph − F lb ) as the net flux from the land surface to the atmosphere (negative for sinks) and S oN = −(F ao − F oa ) as the net flux from the ocean surface to the atmosphere (negative for sinks).After taking S lN = S P lN + S l and S oN = S P oN + S o , where S P lN and S P oN are the prior net CO 2 fluxes to the land and ocean surfaces, respectively, and S l and S o are the residual fluxes to be inverted for the land and ocean surfaces, respectively, Eq. ( 7) can be rewritten as +F lb (δ lb − δ e lb ) + S P oN ε ao + F oa (δ e a − δ a ) .
Equation ( 8) is the theoretical basis for our joint 13 C / 12 C inversion as it links the measured 13 C composition in the atmosphere to the CO 2 fluxes of the land and ocean surfaces.In the implementation of the joint inversion system (Eq.6), a transport matrix is used to link a flux in a particular region to the concentration measured at a particular site.We focus on optimizing the net CO 2 flux using both CO 2 and 13 CO 2 observations rather than optimizing the one-way fluxes, and therefore the discrimination terms to be optimized are moved to the left-hand side of Eq. ( 8) and the disequilibrium terms remain on the right-hand side.Based Eq. ( 8), the regional discrimination D j in Eq. ( 6) is therefore defined as where ε lph,j and ε ao,j are the 13 C fractionation ratio for region j for land and ocean fluxes, respectively.In the joint inversion system, we treat S l and S o as the state variables and D j as predetermined parameters that vary in space (region) and time (monthly).It is therefore a prerequisite to estimate accurately these parameters as well as other isotopic parameters on the right-hand side of Eq. ( 8).
For land regions, BEPS is used to calculate all land variables in Eq. ( 8), including S P lN , F lb , ε lph , R lb , δ lb and δ e lb for each region and month.For ocean regions, ε ao = −2 ‰ and empirical equations developed by Ciais et al. (1995b) are used to calculate F oa and δ e a as functions of sea surface temperature on 1 • × 1 • grids.The 13 CO 2 concentration time series (c m+1 , . ..c m+k ) in Eq. ( 6) in ppm ‰ is the numerical realization of the righthand side of Eq. ( 8) and is computed with the following equation: In Eq. ( 10), C a,i dδ a,i dt can be calculated with observed CO 2 concentration and 13 C composition at two consecutive times, t and t + 1, using the following equation: where C a,i is the mean concentration of CO 2 at each observation station i between t and t + 1, and δ a,i is the 13 C composition at station i, and its derivative with time is taken as its difference between t and t + 1.This derivative represents the δ a growth rate that is the combined outcome of the various isofluxes in Eq. ( 7).The term dt is the sum of 13 δ changes due to fossil fuel and biomass burning, prior land 13 C discrimination flux, land 13 C disequilibrium flux, prior ocean 13 C discrimination flux and ocean 13 C disequilibrium flux, corresponding to the terms in Eq. ( 8). 13 δ k represents the 13 δ value (‰) for each term in Eq. ( 8), and dC k,i dt is the change of concentration (ppm) calculated with the flux of each term in Eq. ( 8) according to the atmospheric transport function M in Eq. ( 6).
The uncertainty of c i as part of the uncertainty matrix R includes the uncertainties of the six terms on the righthand side of Eq. ( 10).The uncertainty for the first term is based on the measurement error (see next Sect.2.1.4)and its global average is 3.08 ppm ‰ month −1 .The uncertainties of terms 2 to 6 are estimated to be 0.95, 3.17, 0.87, 0.12 and 2.69 ppm ‰ month −1 .The total uncertainty for c i is therefore 5.33 ppm ‰ month −1 as a global average, taking as the square root of the sum of the square of the six uncertainties.As an approximation, this total uncertainty is distributed to each station and each month according to the spatial and temporal patterns of uncertainty of the first term.
The inversion system defined by Eq. ( 6) can be implemented in three ways using (1) CO 2 concentration only by excluding the appended matrices for 13 CO 2 , (2) 13 CO 2 data only by using 13 CO 2 -related matrices only and (3) both CO 2 and 13 CO 2 data.Through using the data in these three ways, the information content of 13 CO 2 measurements for CO 2 can be systematically investigated.
In order to investigate the influences of the isotopic discrimination and disequilibrium over land and ocean on the inversion results, we conduct five sets of inversions for the following cases: Case I The spatial variations of all isotopic compositions and the discrimination and disequilibrium fluxes in Eq. ( 8) are considered for both land and ocean.This is the ideal case as the basis to investigate other cases.
Case II The photosynthetic discrimination (ε lph ) over land is taken as a constant of −14.1 ‰, which is the global average obtained by BEPS, and therefore D j = −14.1 ‰.This is a case to ignore regional differences in isotopic discrimination over land.
Case III All isotopic variables are the same as case I, but the land disequilibrium term in Eq. ( 8) is ignored.This is a case to investigate the influence of the land isotopic disequilibrium on the CO 2 flux inversion.
Case IV All isotopic variables are the same as case I, but the ocean disequilibrium term in Eq. ( 8) is ignored.This is a case to investigate the influence of the ocean isotopic disequilibrium on the CO 2 flux inversion.
Case V Both land and ocean disequilibrium terms are ignored, but all other isotopic variables in Eq. ( 8) are same as case I.This is a case to investigate the importance of the total disequilibrium flux in CO 2 flux inversion at the global scale.
Cases III to V are useful not only for evaluating the performance of the joint inversion system but also for assessing the impacts of errors in isotopic disequilibrium estimation on the CO 2 flux inversion.

2.1.4
Covariance matrixes for the CO 2 flux and CO 2 and 13 CO 2 concentration measurements In the joint inversion using both CO 2 and 13 CO 2 measurements, the covariance matrix (Q) for the CO 2 flux remains the same as that in the CO 2 -only inversion (Eq. 3) but the error matrix (R) for concentration measurements is expanded to the dimension of 16 980 × 16 980 to include 60 months of 13 CO 2 observations at 73 stations.Following Deng and Chen (2011), we use an uncertainty of 2.0 Pg C year −1 for the total global land surface CO 2 flux, and this total uncertainty is spatially distributed to the 39 regions according to the annual total net primary productivity (NPP) of these regions simulated by BEPS.For each region, the annual total uncertainty is further distributed to each month according to the simulated seasonal variation in NPP.The global total uncertainty (standard deviation) is spatially and temporally distributed in such a way that the total variance is preserved after the distributions, following the principle of TRANSCOM 3 (Gurney et al., 2003).The uncertainty for the total ocean flux is prescribed as 0.67 Pg C year −1 (Deng and Chen, 2011).
In this way, all the diagonal elements (Q ii ) in the uncertainty matrix Q are determined, while off-diagonal values are assigned to zero, meaning that no flux covariances between regions and months are assumed.The uncertainty of CO 2 measurements in the R matrix is the same as that described in Deng and Chen (2011), following the approach of Peters et al. (2005) and Baker et al. (2006).In this approach, the uncertainty of a monthly CO 2 measurement at a site is estimated as R ii = σ 2 const + GVsd 2 , where constant portion σ const in ppm is assigned according to site category: Antarctic (0.15), oceanic (0.30), land and tower (1.25), mountain (0.90) and aircraft (0.75), while the site-specific variable portion "GVsd" is obtained from the GLOBALVIEW-CO2 2008 database.The 13 CO 2 measurement uncertainty is calculated in a similar way: the variable portion is obtained from the GLOBALVIEW-13CO2 2008 database, while the constant portion is taken as R a σ const in ppm first, where R a is the ratio of 13 CO 2 to CO 2 in the air (∼ 0.011147), and then converted to ‰.The average standard deviation of δ 13 C observations determined in this way for 73 stations is 0.0685 ‰.
2.2 Prior CO 2 and 13 CO 2 flux estimation

CO 2 flux
Terrestrial biosphere fluxes A process-based terrestrial ecosystem model called the BEPS (Chen et al., 1999;Liu et al., 1997) is used in this study to estimate the net terrestrial CO 2 flux and its components including the gross primary productivity (GPP), NPP, heterotrophic respiration (F lb ) and net ecosystem productivity (NEP).GPP is calculated using the Farquhar's leaf-level model (Farquhar et al., 1980) upscaled to the canopy level using a recently refined two-leaf approach (Chen et al., 2012).NPP is taken as 45 % of GPP (Ise et al., 2010) as global biomass data and its components (stem, foliage, root) are lacking for reliable computation of the autotrophic respiration.F lb is calculated as the sum of the decompositional CO 2 release from nine soil carbon pools, namely coarse and dead wood detritus pool, surface structural pool, surface metabolic pool, surface microbial pool, fine-root structural litter pool, fine-root metabolic pool, soil microbial pool, slow carbon pool and passive carbon pool.The sizes of these pools for each cover type in each 1 • grid are estimated using a model spin-up approach based on simulated NPP in 2000 to create a global land sink of 3.73 Pg C year −1 .The total NPP for each 1 • grid is taken as a weighted sum of NPP of seven aggregated land cover types, and the weights are proportional to the areal fractions of the cover types determined using the GLC2000 land cover map at 1 km resolution (Chen et al., 2012).Remotely sensed leaf area index (LAI) (Deng et al., 2006) at 1 km resolution and a clumping index map at 6 km resolution (Chen et al., 2005) and a soil textural map (Webb et al., 1991) are aggregated to 1 • grids for each cover type based on GLC2000 land cover and used as input to BEPS.National Center of Environmental Prediction (NCEP) reanalyzed data (Kalnay et al., 1996;Kanamitsu et al., 2002) are the meteorological drivers for BEPS to simulate hourly carbon fluxes.The output of BEPS used as the prior flux in the inversions is NEP, which does not include carbon emission due to disturbance.

Ocean fluxes
The daily flux of CO 2 across the air-water interface used in this study is constructed based on the results of daily CO 2 fluxes simulated by the OPA-PISCES-T model (Buitenhuis et al., 2006).This model is a global ocean general circulation model (OPA) (Madec et al., 1998) coupled to an ocean biogeochemistry model (PISCES-T) (Aumont et al., 2003;Buitenhuis et al., 2006).PISCES-T represents the full cycles of C, O 2 , P, Si, total alkalinity and a simplified Fe cycle.It also includes a representation of two phytoplankton, two zooplankton and three types of dead organic particles of different sinking rates.OPA-PISCES-T is forced by daily wind stress and heat and water fluxes from the NCEP reanalyzed data (Kalnay et al., 1996;Kanamitsu et al., 2002).Hourly S o ( 13 C) is calculated with gridded optimum interpolation sea surface temperature from the NOAA National Climate Data Center (Reynolds and Smith, 1994;Reynolds et al., 2002).

Fossil-fuel emissions
The fossil-fuel emission field (2000)(2001)(2002)(2003)(2004) used in this study (http://carbontraacker.noaa.gov) is constructed based on (1) the global, regional and national fossil-fuel CO 2 emission inventory from 1871 to 2006 (CDIAC) (Marland et al., 2009) and (2) the EDGAR 4 database for the global annual CO 2 emission on a 1 • grid (Olivier et al., 2005).The 13 CO 2 flux from fossil-fuel consumption is calculated from CO 2 emissions of different fuel types multiplied by their respective 13 C / 12 C ratios with consideration of their latitudinal distributions based on Andres et al. (2000).

Fire emissions
CO 2 emissions due to vegetation fires are an important part of the carbon cycle (van der Werf et al., 2006).Each year, vegetation fires emitted around or more than 2 Pg C of CO 2 into the atmosphere, mostly in the tropics.The fire emission field used in this study is based on the Global Emissions Fire Database version 2 (GFEDv2) (Randerson et al., 2007;van der Werf et al., 2006).

13 CO 2 flux
Based on the initial work of Chen et al. (2006), BEPS is further developed to include a capacity to compute the global distribution of the terrestrial 13 CO 2 flux.Following the principle of multi-stage 13 C fractionation in the pathway through leaf boundary layer, stomates, messophyll and chloroplast initially proposed by Farquhar et al. (1984Farquhar et al. ( , 1989) ) and implemented globally by Suits et al. (2005), we developed a module in BEPS for computing the total photosynthetic fractionation and the resultant 13 CO 2 flux.Specifically, the photosynthetic discrimination for C3 plants ( PC3 ) is calculated from where b , s , diss , aq and f are the rates of discrimination against 13 CO 2 through leaf boundary layer, stomates, dissolution in mesophyll water, transport in aqueous phase and fixation in chloroplast, respectively, and are assigned values of 2.9, 4.4, 1.1, 0.7 and 28.2 ‰ (Suits et al., 2005).
A is the photosynthetic rate in mol m −2 s −1 and p equals 0.022624 T a /(273.16P ) with the dimension of m 3 mol −1 , where T a is air temperature in • K and P is the standard air pressure at 1.013 bar.C a and C c are the CO 2 concentrations in mol mol −1 in the free air and leaf chloroplast, respectively.For C4 plants, the photosynthetic discrimination ( PC4 ) is taken as a constant of 4.4 ‰ (Suits et al., 2005).
The leaf boundary layer (g b ) is calculated with the following equation: where α is the diffusivity of CO 2 in dry air in m 2 s −1 calculated as 10 −6 (0.129 + 0.007 T a ) and T a is the air temperature in • C; l is the leaf characteristic dimension in metres, taken as a constant of 0.1 m; and N is the Nusselt number equal to (u d l/υ) 0.5 , where u d is the wind speed in m s −1 at the vegetation displacement height (80 % of the average vegetation height) and υ is the kinematic viscosity of dry air in m 2 s −1 calculated as 10 −6 (0.133 + 0.007 T a ).u d is derived from the wind speed above the canopy based on LAI and vegetation height assigned according to plant functional type (Table 1).
As part of the GPP calculation, the stomatal conductance (g s ) computed separately for sunlit and shaded leaves using the Ball-Berry equation (Ball, 1988), where f w is a scaling factor depending on soil moisture and texture (Chen et al., 2012); h s is the air humidity at the leaf surface; C s is the CO 2 concentration at the leaf surface; p is the same as in Eq. ( 12); and m and b are the slope and intercept in this linear relationship, and they are assigned values according to plant function type (Table 1) (Chen et al., 2012).The mesophyll conductance g m is calculated based on the method of Harley et al. (1992): where A is the photosynthetic CO 2 assimilation rate; C i is partial pressure of CO 2 in the air spaces inside leaves; R d www.geosci-model-dev.net/10/1131/2017/ is the respiration rate occurring during the day not related to photorespiration; is the CO 2 compensation point in the absence of R d ; and J is the rate of photosynthetic electron transport.These parameters are the same as those used in computing the CO 2 flux.
Our methods of computing stomatal and mesophyll conductances differ from previous studies (Suits et al., 2005;Scholze et al., 2008;Rayner et al., 2008) in the following ways: (1) these conductances are calculated separately for sunlit and shaded leaves because BEPS is a two-leaf model, in which the total GPP of a canopy is taken as the sum of sunlit and shaded leaf GPP, and (2) the mesophyll conductance mechanistically depends on a set of parameters rather than being treated as a constant or to be proportional to the stomatal conductance.Since it has been demonstrated that sunlit and shaded leaf separation is essential for accurate modelling of canopy-level photosynthesis (Chen et al., 1999;Sprintsin et al., 2011), it is expected that this separation is also essential for 13 CO 2 flux estimation.We found that the use of Harley's method for computing the mesophyll conductance makes the calculated 13 C photosynthetic fractionation stable for its global application, while the simpler method of treating the mesophyll conductance in proportion with the stomatal conductance often incurs abnormally large or small values of 13 C photosynthetic fractionation.
The photosynthetic 13 CO 2 flux is in disequilibrium with the respiratory 13 CO 2 flux because of the change in atmospheric 13 CO 2 concentration since the preindustrial time (Ciais et al., 1995b;Fung et al., 1997).The heterotrophic respiratory flux from the decomposition of organic matter of different ages carries the memory of the past atmospheric 13 CO 2 concentration, while the photosynthetic 13 CO 2 flux is affected by the current atmospheric 13 CO 2 concentration.The isotopic composition of each of the nine soil carbon pools (δ 13 C soil,i ) is estimated with following formula: where δ 13 C a is the isotopic composition of carbon in atmospheric CO 2 in the past as determined by the ice-core record (Francey et al., 1999), ε lph is the annual mean of photosynthetic discrimination in 2003 and τ i is the age of carbon pool i (Table 2) (Ju and Chen, 2005).In the calculation of the mean age of a carbon pool, we have considered the ages of various carbon pools at the time of entering the pool (Potter et al., 1993) so that the mean age is considerably larger than the turnover time determined by the decomposition rate (Fung et al., 1997).The mean δ 13 C soil is taken as the flux-weighted δ 13 C soil,i for the nine carbon pools.The results of δ 13 C soil for the globe are shown in Fig. 5.The 13 C composition of the biosphere δ lb in Eq. ( 8) is taken as the mean δ 13 C soil , while the biospheric 13 C composition δ e lb in equilibrium with the current atmosphere is taken as δ a − ε lph .
The accuracy of the BEPS model in simulating atmospheric 13 CO 2 concentration was previously tested (Chen et al., 2006;Chen and Chen, 2007) against measurements over a boreal forest at Fraserdale, Ontario, Canada (49 • 52 29.9 N, 81 • 34 12.3 W).Flask measurements of δ 13 C a were made 40 times in both daytime and nighttime on a tower at a height of 20 m during a 3-day campaign on 21-23 July 1999.BEPS simulated these measurements with RMSE = 0.34 ‰ and r 2 = 0.76.

Transport modelling
A transport-only version of the atmospheric chemistry and transport model TM5 (Krol et al., 2003(Krol et al., , 2005) ) is used for CO 2 and 13 CO 2 transport modelling to produce a fully linear operator on these fluxes.The spatial resolution of TM5 is 6 • × 4 • for the globe and 3 • × 2 • for North America, and the atmosphere is divided vertically into 25 layers with 5 layers in the planetary boundary layer.Tracer transport (advection, vertical diffusion, cloud convection) in TM5 is driven by offline meteorological fields taken from the European Centre for Medium Range Weather Forecast (ECMWF) model.All physical parameterizations in TM5 are kept the same as the ECMWF formulation to achieve compatibility between them.The four background fluxes from terrestrial ecosystems, oceans, fossil-fuel burning and biomass burning are individually inputted to TM5 to calculate the contributions of these fluxes to the atmospheric CO 2 and 13 CO 2 concentrations.Since the main purpose of this study is to develop a joint inversion system, only one transport model is used, the transport matrix M is assumed to be free of errors.

CO 2 and 13 CO 2 data sets
Monthly CO 2 and 13 CO 2 concentration data from 2000 to 2004 are compiled from the GLOBALVIEW CO 2 and 13 CO 2 database.Though the GLOBALVIEW database consists of both extrapolated and interpolated data that were created based on the technique devised by Masarie and Tans (1995), we selected the synchronized and smoothed values of actual observations to compile our concentrations data sets.Only direct measurements of CO 2 from the GLOBALVIEW data set are used in our inversion after using a time-frequency weighting scheme (Deng and Chen, 2011).There are 5431 monthly data from 209 sites for 42 months used for CO 2 (5431 out of 8778, i.e. 209 × 42), and 3066 monthly data from 73 sites for 13 CO 2 (i.e.73 × 42 monthly data).Since the number of 13 CO 2 observation sites is much smaller than that of CO 2 sites, all monthly data at 73 sites are used for 13 CO 2 , and the missing 13 CO 2 data are filled with the reference data provided in the same GLOBALVIEW data set.The filled data may have introduced an additional error to the data set as shown in Fig. 15b.
To minimize the non-linear aggregation effects of the large regions (Pickett-Heaps, 2007), the contributions of the four background fluxes are subtracted from the above monthly concentrations.So the matrix c in Eqs. ( 3) and ( 4) is expressed as where c obs is the monthly CO 2 and 13 CO 2 concentrations obtained from GLOBALVIEW, and c ff , c bio , c ocn and c fire are simulated contributions of CO 2 and 13 CO 2 concentrations from the terrestrial biosphere, ocean, fossil-fuel and fire fluxes, respectively.

Prior CO 2 and 13 CO 2 fluxes
Terrestrial ecosystem models integrate many sources of information, including vegetation structure, soil and meteorology, to estimate carbon exchange of the land surface with the atmosphere.Prior CO 2 and 13 CO 2 fluxes produced by a model can therefore provide indispensable constraints to the otherwise ill-posed inversion based on CO 2 and 13 CO 2 concentration observations alone.Depending on the assigned relative magnitudes of the error matrixes of these observations and these prior fluxes (i.e.R and Q in Eq. 3), these prior fluxes can have equal or even dominant importance to these observations in the inversion results.We have therefore paid a great attention in modelling these prior fluxes, in order to minimize the total inversion errors.Figure 2a shows an example of the global terrestrial GPP distribution in 2003 modelled by BEPS.The total GPP in this year is 132 ± 22 Pg C year −1 (Chen et al., 2012).This value is larger than some of the recent estimates, such as 123 Pg C year −1 by Beer et al. (2010), mostly because the LAI values used as input to BEPS are generally larger than those of the MODIS product (Garrigues et al., 2008).Our LAI values are larger because we used a global clumping index map derived from a multi-angle satellite sensor POLDER (Chen et al., 2005).Clumping increases shaded leaves, which contributed about 35 % to the total GPP globally.Without considering this clumping effect, the shaded leaf area is underestimated resulting in an underestimation of the global GPP by 9 % (Chen et al., 2012).As the spatial distribution of clumping is not uniform (boreal and tropical forests are most clumped and crops and grasses are least clumped), this refinement in the GPP spatial distribution would have some effects on the inversion results between regions.
The NEP, which is the difference between GPP and ecosystem respiration modelled by BEPS, is shown in  smaller (globally 2 Pg C year −1 by BEPS) because a model spin-up approach is used to estimate the soil carbon pool sizes based on a dynamic equilibrium assumption.Under this assumption, the annual heterotrophic respiration (F lb ) equals annual NPP during the preindustrial period, and the soil carbon pool sizes are derived from F lb by solving a set of differential equations describing the decomposition and interactions among the pools (Govind et al., 2011).In this way, F lb is forced to depend on NPP and the systematic biases in GPP are not carried into NEP estimation.NEP is non-zero after the preindustrial period because of the changes in climate and atmospheric composition (CO 2 and nitrogen) as well as disturbance.In our regional modelling, both disturbance and non-disturbance effects are considered for Canada (Chen et al., 2003) and USA (Zhang et al., 2012) forests.However, in our global model spin-up from 1901 (taken as the end of preindustrial period) to 2000, only the non-disturbance effects are considered because of lack of spatially explicit disturbance data outside of North America, while carbon emission due to fire disturbance in the study period from 2000 to 2004 is considered separately using the GFED data set (Randerson et al., 2007;van der Werf et al., 2006).The prior net CO 2 fluxes for the globe for the years 2002-2004 are given in Table 3 with inversion results with and without the 13 C constraint.
Table 3. Inverted fluxes (Pg C year −1 ), averaged for 2002-2004, for land and ocean regions with (CO 2 + 13 CO 2 ) and without (CO 2 -only) 13 C constraint.The negative sign denotes the flux from the atmosphere to the surface (sink).Various treatments are made to 13 C discrimination and disequilibrium fluxes represented by the following cases.Case I: full consideration of the regional differences in discrimination and disequilibrium; case II: same as case I, but the annual photosynthetic discrimination ratio is set at a constant of −14.1 ‰, although it is monthly variation pattern as modelled by BEPS is retained; case III: same as case I, but the disequilibrium flux over land is ignored; case IV: same as case I, but the disequilibrium flux over ocean is ignored; case V: same as case I, but the disequilibrium flux over both land and ocean is ignored.
Inverted  The global distribution of the total photosynthetic discrimination (δ 13 C pt = δ 13 C a − ) modelled by BEPS is shown in Fig. 3. Forests, such as those in North America, Russia, Europe, Amazon, central Africa, central China and southeast Asia, generally have high photosynthetic discrimination rates (> 16 ‰), whereas grassland and cropland (in particular C4 grasses and crops) have low discrimination rates.Also shown in Fig. 3 is the ocean diffusive discrimination against 13 CO 2 .The discrimination over ocean is much smaller than that over land.This difference between land and ocean discrimination may be considered as the largest signal of 13 CO 2 observations on the global carbon cycle (Tans et al., 1990;Rayner et al., 2008) and is considered in our inversion using different 13 CO 2 discrimination rates for ocean and land regions (see Eq. 6).
To estimate the disequilibrium between photosynthetic and respiratory discrimination against 13 CO 2 , the global dis-tribution of the mean soil carbon age is computed after weighting the ages of the nine soil carbon pools against their fluxes due to decomposition (Fig. 4).Forests at high latitudes have the soil carbon age of about 40-60 years, while the tropical forests have much lower values in the range from 10 to 30 years.This latitudinal distribution pattern is mostly determined by soil temperature.In low latitudes, high temperature induces fast turnovers of detritus and fast soil carbon pools, whereas at high latitudes, low temperature maintains relatively large fractions of slow and passive soil carbon pools.Cropland and grassland also have larger fractions of fast and detritus carbon pools than forest cover types and therefore have younger soil carbon on average.This spatial distribution of soil carbon age has a strong influence on the total respiratory discrimination against 13 C (δ 13 C r ) calculated by BEPS (Fig. 5).Respiration from older carbon at high latitudes carries the memory of the older atmosphere with   less 13 CO 2 concentration and hence has lower discrimination rates (larger δ 13 C r or smaller absolute value).However, respiration would mostly depend on the photosynthetic discrimination rates as soil organic matter originates from photosynthetic production.As a result, forested areas have higher respiratory discrimination rates (lower δ 13 C r or larger absolute value).Most of the high values of δ 13 C r in Fig. 5 are associated with large fractions of C4 plants in the grid, such as the corn belt in the USA, cropland in northeastern China, southern border of Sahara desert and southeastern South America.
The global distribution of the disequilibrium between photosynthetic and respiratory discrimination, taken as the difference between Fig. 3 and 5, is shown in Fig. 6.The disequilibrium is the largest at high-latitude boreal forests in North America and Eurasia because their soil carbon is the oldest, as shown in Fig. 4. The spatial distribution pattern of the disequilibrium is similar to those of Ciais et al. (1995b) and Fung et al. (1997) but the magnitude is larger because the date of our result in 2000 are more recent than these two previous studies.As the time lapses, the atmosphere is getting lighter in terms of the isotopic composition of CO 2 resulting from the increased air-borne CO 2 from fossil-fuel consumption.Also shown in Fig. 6 is the disequilibrium over the ocean estimated using the method of Ciais et al. (1995b).This ocean disequilibrium has a large latitudinal gradient because of the gradients in sea surface temperature gradient and the fluxes of CO 2 and 13 CO 2 .The spatial distribution in the disequilibrium and the differences in disequilibrium between ocean and land may be considered to be the secondary signal of 13 CO 2 observations on the global carbon cycle.The effects of these disequilibria on the carbon flux are considered in our inversion through pre-subtracting their contributions to the measured 13 CO 2 composition in Eq. ( 10).

Inverse modelling results
Although the inversions were made for the 2000-2004 period, the results of the first 2 years are not included in the analysis because they are affected by the assumption of uniform global distributions of CO 2 and 13 CO 2 concentrations at the start of our transport modelling using TM5.An 18-24month period is usually considered to be necessary for the simulated distributions to reach realistic states with reasonably accurate prior surface fluxes from ocean and land and atmospheric transport simulations (Rödenbeck et al., 2003;Deng and Chen, 2011).The following results are therefore summarized as the average for the 2002-2004 period.

Partition between ocean and land sinks with and without 13 CO 2 constraint
To investigate the usefulness of 13 CO 2 observations in partitioning between ocean and land sinks, we conducted inversions with and without 13 CO 2 constraint as expressed in Eq. ( 6), i.e. with and without the 13 C-related expansions of the matrixes.The CO 2 -only inversion increases the land sink from the prior of 2.61-3.40Pg C year −1 while decreasing the ocean sink from the prior of 2.13-1.48Pg C year −1 (Table 3).These results are similar to those of Deng and Chen (2011).
The results from the joint inversion are considerably differ-ent; the posterior sinks for land and ocean become 2.53 and 2.36 Pg C year −1 (Table 3), respectively, suggesting that the use of 13 CO 2 observations in the inversion considerably influenced the partition between land and ocean fluxes.The ratio between land and ocean sinks is 1.07.The joint inversion system developed in this study may be regarded as a different form of double deconvolution.Using the double deconvolution method with the global average disequilibrium coefficients of 0.49 and 0.78 ‰ and the disequilibrium fluxes of 26.8 and 66 Pg C year −1 ‰ for land and ocean derived in this study (Table 4), respectively, we also calculated the land and ocean sinks to be 2.90 and 2.36 Pg C year −1 , respectively.The ratio between land and ocean sinks is 1.23, which is close to the value of 1.07 derived from the joint inversion system, indicating that the joint inversion can effectively perform double deconvolution.Our joint inversion system differs from previous double deconvolution systems (Siegenthaler and Oeschger, 1987;Keeling et al., 1989a;Francey et al., 1995;Randerson et al., 2002) in the following ways: (1) the estimation of CO 2 fluxes for the land and ocean is additionally constrained by the prior fluxes for the land and ocean rather than being entirely dependent on measured CO 2 concentration and 13 CO 2 composition, and (2) the spatiotemporal variations in all parameters associated with isotopic discrimination and disequilibrium are considered in the estimation of the CO 2 flux using a mechanistic biospheric model rather than global average values or simple models based on covariates.These differences in methodology as well as the differences in the mean disequilibrium fluxes may explain why the ocean and land sinks from the joint inversion system differ from the various double deconvolution results.
3 * F f is the carbon emission from fossil fuel and biomass burning, 6.9 and 2.1 Pg C year −1 , respectively, and δ f is weighted average 13 C composition for fossil fuel and biomass burning, being 25.27 ‰, and δ a = −8.0‰.
The impacts of 13 CO 2 data on the joint inversion can also be evaluated from the view point of global 13 CO 2 mass budget.Table 5 shows the budgets and its components for the prior, double deconvolution, CO 2 -only inversion and joint inversion cases.In these cases, the isofluxes due to fossil-fuel emission, land and ocean disequilibrium and atmospheric storage change are the same, and only those due to discrimination over land and ocean are adjusted.The prior case shows a global imbalance of −5.0 Pg C year −1 ‰, indicating that either the prior land or ocean fluxes or both are inconsistent with 13 CO 2 measurements.Through double deconvolution, this imbalance is greatly reduced to 0.8 Pg C year −1 ‰, mostly by an increase in the discrimination flux over land because of its large discrimination rate.The CO 2 -only inversion increases the land discrimination flux while decreasing the ocean discrimination flux, resulting in no improvement in the global isotopic balance.The joint inversion optimized both ocean and land fluxes in the direction consistent with 13 CO 2 measurements, reducing the imbalance considerably to 1.8 Pg C year −1 ‰.These cases illustrate clearly that the global isotopic mass balance is very sensitive to the partition between ocean and land fluxes because of the large differ- ence in the discrimination rate between land and ocean.In this analysis, the disequilibrium fluxes are not adjusted, but the influences of the uncertainties in these fluxes on the inversion results are analyzed in Sect.3.2.4.
Existing estimates for the ocean sink for anthropogenic CO 2 in the 2000s varies from 1.94 to 2.6 Pg C year −1 (Wanninkhof et al., 2013;Landschützer et al., 2014;Majkut et al., 2014;DeVries, 2014).The average ocean sink for the 2002-2004 period summarized by the Global Carbon Project (GCP) (Le Quéré et al., 2013) is 2.4 Pg C year −1 , while the land sink in the same period is 2.7 Pg C year −1 as the residual of the global carbon budget after including the emission due to land use change as a source of carbon.Although the prior estimates of these sinks in our inversions are similar to these values, our CO 2 -only inversion considerably increases the land sink and decreases the ocean sink.The addition of 13 CO 2 measurements to the inversion significantly decreases the land sink and increases the ocean sink, pulling the inversion results in the direction to agree with these existing estimates (Fig. 7).This may indicate that the use 13 CO 2 measurements in the joint inversion has improved the CO 2 estimation.In this comparison, we have not considered the unknown small amount (0.1-0.3 Pg C year −1 ) of lateral carbon transport in rivers from land to ocean.This amount is included in some of the estimates of the ocean sink used by GCP, and therefore should be subtracted from the ocean sink and added to the land sink by GCP in order to compare with our atmospheric inversion results.

Influence of 13 CO 2 constraint on the spatial distribution of the inverted carbon flux
The 13 CO 2 constraint modified not only the partition between ocean and land fluxes but also their spatial distribution patterns.Figure 8 shows the result of the CO 2 -only inversion (i.e.without the 13 CO 2 constraint), as the net carbon flux over land and ocean averaged for the period of 2002-2004.
Figure 9 shows the difference between inversions with and without the 13 CO 2 constraint, i.e. the result of CO 2 + 13 CO 2 inversion minus that of CO 2 -only inversion.The general patterns of the inverted carbon flux are similar between these two inversions because these inversions depend primarily on the CO 2 concentration, the prior flux, the error matrixes of the prior flux, and concentration observations.However, there are several large or notable differences: (1) the Amazon region (region 31) is changed from a carbon source to a sink (Fig. 10; note: a reduction in sources is shown as a negative value); (2) the carbon sink in the tropical Asia (region 37) is noticeably reduced (by about 10-20 g C m −2 year −1 from a sink magnitude of about 80-100 g C m −2 year −1 ); (3) the sink in Asia (region 36) decreases pronouncedly by about 10-20 g C m −2 year −1 , while the sinks in Russia (region 35) and Europe (region 39) are also reduced by some extents (about 5-20 g C m −2 year −1 ); (4) most small regions in the southern part of North America show increases in sinks, but those in the northern part (Canada and Alaska) show increases in sources (see also Fig. 11); the overall sink in North America decreases from 0.67 to 0.54 Pg C year −1 (Fig. 10); and (5) most ocean regions at mid-latitudes have small gains in sink.
It is of particular importance to note that the 13 CO 2 constraint changed the Amazon region from a carbon source of 0.43 ± 0.46 Pg C year −1 to a carbon sink of 0.08 ± 0.38 Pg C year −1 with a notable reduction in the posterior uncertainty, which is higher than uncertainty reductions in most other regions (Fig. 10).This change is likely caused by the relatively large addition of information from 13 CO 2 in this tropical region where CO 2 observations are sparse, causing large uncertainties in the inverted flux in this region in the CO 2 -only inversion.Potter et al. (2009) simulated the NEP of the Amazon region using the Carnegie-Ames-Stanford Approach driven by remote sensing inputs and found that the NEP for the region was slightly negative (−0.07 Pg C year −1 ) over the 2000-2004 period.Davidson et al. (2012) summarized from various inventory-based studies that mature forests in the region were accumulating carbon at a rate of 0.29-0.57Pg C year −1 over the decade before 2005, meaning that NEP is positive.Since the fire emission is estimated to be 0.50 Pg C year −1 (Richey et al., 2002), the Amazon region would be either a net source of carbon or about carbon neutral.Since spatially explicit fire emission is considered together with fossil-fuel emission as a source in our study, the inverted carbon flux corresponds to −NEP, and therefore the result from our joint inversion is in broad agreement with the results of Potter et al. (2009) and Davidson et al. (2012).Without the 13 CO 2 constraint, our inversion result shows an unreasonably large source of carbon in the Amazon region.

Influence of the spatial distribution of photosynthetic discrimination on the inverted carbon flux
The joint inversion results shown in Figs.9-11 are from case I with the best estimates of the 13 C discrimination and disequilibrium fluxes and therefore represent a baseline study to which other cases are compared for the purpose of investigating the importance of accurate consideration of the spatial distributions of isotopic discrimination and disequilibrium over land and ocean.Case II is designed to investigate the importance of considering the spatial distribution of the photosynthetic isotopic discrimination over land for inverting the CO 2 flux by fixing the discrimination at a constant over land.Figure 12a shows the spatial distribution of the difference in the total isotopic discrimination, i.e.D j = ε lph,j , among 39 land regions between case I and case II, calculated as case I minus case II.Regions with positive differences in D j are shown with positive differences in the inverted CO 2 flux (Fig. 12b), meaning larger sinks (negative values) in case II, and vice versa.This is because a smaller discrimination rate (smaller than −14.1 ‰) means a larger CO 2 flux from the atmosphere to the surface (more negative value) for the same change in 13 CO 2 concentration in the  atmosphere.Under the same condition, a larger discrimination induces a smaller sink (less negative).The absolute regional differences between case I and case II are considerable (Fig. 12b), e.g. up to 18 g C m −2 year −1 , showing increases in sinks in Africa, Asia and Australia and decreases in sinks in Amazon, Europe, Russia and most of the small regions in North America.However, the total global sink values of case II after ignoring the spatial distribution of the disequilibrium rate over land change very little those of case I (Table 3): from 2.53 ± 0.93 to 2.49 ± 0.95 Pg C year −1 for land and from 2.36 ± 0.49 to 2.35 ± 0.48 Pg C year −1 for ocean.This is because the global mean discrimination rates are the same between these two cases.

Influence of the uncertainties in disequilibrium fluxes on the inverted carbon flux
The average disequilibrium coefficients and fluxes for land and ocean derived in this study are comparable to published results (Table 4), although the estimates of the disequilibrium flux over ocean in previous studies vary in a large range.The uncertainty in the estimated land and ocean disequilibrium fluxes mainly arises from two sources: the estimated disequilibrium coefficient and one-way CO 2 flux from the surface.Mathematically, the total uncertainty in the disequilibrium flux, denoted as For land, the first source depends  Case III, case IV and case V are conducted to investigate the relative importance of the disequilibrium fluxes over land and ocean (Table 3) in the CO 2 flux inversion.In case III, where the disequilibrium over land is ignored while other settings remain the same as case I, the land sink increases by 1.05 Pg C year −1 , whereas the ocean sink decreases by 0.08 Pg C year −1 in comparison with case I.When the disequilibrium over ocean is ignored instead (case IV), the land sink increases by 0.13 Pg C year −1 , while the ocean sink increases by 2.08 Pg C year −1 , in comparison with case I.When the disequilibria over both land and ocean are ignored, the land sink increases by 1.18 Pg C year −1 , while the ocean sink increases by 1.96 Pg C year −1 , in comparison with case I. Results from these case studies suggest that in the joint inversion using both CO 2 and 13 CO 2 measurements, the inverted CO 2 flux can be significantly influenced by the disequilibrium fluxes of land and ocean.The carbon sinks over land and ocean increase when these disequilibrium fluxes are ignored because the photosynthetic and diffusive sources of 13 CO 2 have to increase to make up for the shortfall due to ig- noring the disequilibrium sources.These pronounced influences of the disequilibrium on the CO 2 sink inversion suggest that 13 CO 2 data contain strong signals for the global carbon cycle.In the joint inversion, these data can have the power to distort the global CO 2 mass balance if the 13 CO 2 mass budget (Eq.8) is not properly simulated.The influence of 13 CO 2 on the joint inversion depends only weakly on the estimated uncertainty in the 13 CO 2 data.We found that if the uncertainty is reduced by half, the sum of the land and ocean sink deviates from the CO 2 -only case by 2-6 % for all scenarios, suggesting that the mean disequilibrium fluxes play the dominant roles in the joint inversion.
The impacts of these disequilibrium fluxes on the inverted CO 2 flux determined in case III, case IV and case V are similar to previous results using the double deconvolution tech-nique (Tans et al., 1993;Ciais et al., 1995b;Randerson et al., 2002).However, the influences of these disequilibrium fluxes on the joint inversion could possibly be compromised due to the small number of 13 C observation sites relative to the number of CO 2 observation sites used in the joint inversion.The number of linear equations for CO 2 concentration in our joint inversion system (Eq.6) greatly exceeds the number for 13 C composition, with a potential of dampening the impact of 13 C data on the inverted results.To investigate the possibility of this dampening effect, we conducted a set of inversions using 13 C data alone (Table 6) and found that the impacts of the disequilibrium fluxes on the inversion results are similar to those of the joint inversion.In case V shown in Table 6, for example, ignoring the disequilibrium fluxes causes the land sink to increase by 1.06 Pg C year −1 and ocean sink to in-Table 6. Inverted fluxes (Pg C year −1 ), averaged for 2002-2004, for land and ocean regions using 13 C data only.The negative sign denotes the flux from the atmosphere to the surface (sink).Various treatments are made to 13 C discrimination and disequilibrium fluxes represented by the cases outlined in Table 3 crease by 2.37 Pg C year −1 , resulting in a total increase of 3.43 Pg C year −1 , which is similar to the total difference of 3.14 Pg C year −1 produced by the joint inversion.These similar results suggest that 13 C data used in the way described by Eqs. ( 6)-( 8) have played the expected role in the joint inversion.By comparing results shown in Tables 3 and 6, it is also encouraging to see that inversions using 13 C data alone can produce reasonable results for the CO 2 flux, although we believe that the joint inversion results shown in Table 3 are more reliable.Our findings on the usefulness of the small 13 CO 2 data set somewhat confirms the claim of Enting et al. (1993Enting et al. ( , 1995) ) that the temporal trend in 13 CO 2 concentration is the major signal constraining the partition between ocean and land sinks.
According to the difference of the inverted flux between case III to case I, the uncertainty of 8.0 Pg C year −1 ‰ in the land disequilibrium flux would cause an uncertainty of 0.47 Pg C year −1 in the land flux.According to the comparison between case IV to case I, the uncertainty of 12.7 Pg C year −1 ‰ in the ocean disequilibrium flux would cause an uncertainty of 0.54 Pg C year −1 in the ocean flux.These uncertainties in the land and ocean fluxes are 17 and 24 % of the jointly inverted fluxes for land and ocean (case I in Table 3), respectively.The impact of the uncertainty in the disequilibrium flux over land is only slightly smaller than the posterior uncertainty of the inverted land flux, but the impact over ocean is larger than the posterior uncertainty.

Discussion
After the CO 2 fluxes are optimized through the inversions, the posterior CO 2 concentration at all stations in each month can be calculated from Eq. (2), and similarly the posterior 13 CO 2 composition can also be calculated from Eq. ( 10) by replacing the prior discrimination fluxes with posterior discrimination fluxes.One way to evaluate the effectiveness of the joint inversion is to examine the improvement in the posterior CO 2 and 13 CO 2 concentrations against measurements.Figure 13 shows concentrations for 10 randomly selected stations from different regions, which are indicated in Fig. 1.The CO 2 and 13 CO 2 concentrations produced using the prior fluxes considerably deviate from observations at all stations.
The posterior CO 2 concentration from the CO 2 -only inversion shows great improvements over the prior concentration in comparison with observations.The posterior CO 2 concentration from the joint inversion does not differ significantly from that of the CO 2 -only inversion.At some stations, the joint inversion produces slightly lower root mean square differences (RMSD) against observations, but in some stations the opposite is true, as indicated by the RMSD values shown in the header of each plot.It is expected that in some stations, the posterior CO 2 concentration in the joint inversion can be slightly worsened because of the influence of 13 CO 2 .The posterior 13 CO 2 concentration is pronouncedly improved over the prior in comparison with observations and almost mimics the observed magnitudes and temporal variations, indicating that the joint inversion system can forcefully adjust CO 2 fluxes to match with 13 CO 2 observation through the prescribed discrimination rates.The posterior CO 2 concentrations for either CO 2 -only or joint inversion show larger seasonal amplitudes than observations at Northern Hemisphere stations, although the means are about the same as observations.This suggests that both carbon uptake during the growing season and ecosystem respiration in the nongrowing season might have been overestimated, even though the annual net carbon flux may be unbiased.Further work is needed to constrain the large photosynthetic and respiratory fluxes separately rather than the net flux only.
In order to provide a comprehensive evaluation, the posterior CO 2 and 13 CO 2 concentrations at all stations are shown in Figs. 14 and 15 against observations.In Fig. 14, we see pronounced improvements in the posterior concentrations from both the CO 2 -only and joint inversions over the prior case.However, the improvements of these two inversions are similar (the joint inversion has a smaller intercept and a slope closer to one, but the CO 2 -only inversion has a slightly larger r 2 value).This is in agreement with the cases shown for the individual stations: some stations are improved and some worsened by the use of 13 CO 2 data, manifesting the force of this additional data on the inversion.In Figure 15, the posterior 13 CO 2 concentration from the joint inversion is shown to be greatly improved from the prior case.In the joint inversion, the increase of the posterior land and ocean sinks over the prior sinks that remove CO 2 from the atmosphere Figure 13.Left panel: comparison of CO 2 concentrations calculated using the prior flux (solid red) and from CO 2 -only inversion (dashed purple) and joint inversion (dashed green) against observations (blue) at 10 randomly selected stations from different regions.The header of each plot indicates the station ID and the root mean square difference (RMSD) for the prior, joint and CO 2 -only inversions against observations.Right panel: comparison of 13 CO 2 composition from the prior (solid red) and joint inversion (dashed green) against observations (blue).The header of each plot indicates the station ID and RMSD of the prior and the joint inversion against observations.logically corrects for the positive bias in the CO 2 concentration produced using the prior fluxes (Fig. 14a).The posterior concentration correlation with observation is stronger for 13 CO 2 than for CO 2 , indicating that isofluxes are effectively optimized in the joint inversion according to 13 CO 2 data.However, some points in Fig. 15b scatter greatly from the 1 : 1 line, and these points are mostly likely the missing data filled with the reference data (Sect.2.4).As other error sources cannot be excluded, these data are retained in our inversion.
After adding 13 CO 2 data to the inversion system, the uncertainty in the inverted CO 2 flux increased from 0.84 to 0.93 Pg C year −1 for land and from 0.40 to 0.49 Pg C year −1 for ocean (Table 3, difference between the CO 2 -only case and case I); i.e. 11 and 23 % increases in uncertainty for land and ocean, respectively.The relative error in preprocessed 13 CO 2 measurements used in the joint inversion is considerably larger than that in CO 2 measurements, causing these increases in the uncertainty of jointly inverted CO 2 fluxes from the CO 2 -only case.The 13 CO 2 measurements were preprocessed before the inversion as the remaining concentration after removing the contributions of fossil-fuel emission and prior land and ocean discrimination and disequilibrium fluxes (Eq.10), and therefore they contain uncertainties from these contributions in addition to measurement uncertainties.Errors in modelling the spatial and temporal variations of the 13 CO 2 flux stem from many sources including errors in modelling the discrimination, which is affected by the fractionation of the 13 CO 2 flow through leaf boundary layer, stomata, mesophyll, etc., and the disequilibrium, which depends on the sizes of nine soil carbon pools and their ages.Although the ocean 13 CO 2 discrimination is small, its disequilibrium has a strong latitudinal gradient, which is approximately calculated using the mean monthly temperature.The error in the calculated ocean disequilibrium coefficient is estimated to be ±1.2 ‰ for the monthly values at a given location and ±0.12 ‰ for the global annual total.Because of these errors, we estimate that the relative uncertainty in the prior 13 CO 2 fluxes is similar to that of the prior CO 2 flux over both land and ocean.

Conclusions
The usefulness of atmospheric 13 CO 2 measurements at 73 stations for global carbon cycle estimation is explored through their use as an additional constraint on an atmospheric inversion of the surface carbon flux using CO 2 observations.The following conclusions are drawn from this study: 1.This 13 C constraint on the joint inversion considerably alters the partition between land and ocean sinks obtained from CO 2 -only inversion, decreasing the land sink from 3.40 ± 0.84 to 2.53 ± 0.93 Pg C year −1 , while increasing the ocean sink from 1.48 ± 0.40 to 2.36 ± 0.49 Pg C year −1 for the 2002-2004 period.Over land, this alteration induces the largest sink increases in the Amazon region and the largest source increases in southern Africa, and Asia, where CO 2 observations are sparse and therefore the additional signal from 13 CO 2 data becomes most important.Over the ocean, sink increases are found broadly at mid-and high latitudes in both hemispheres.
2. The spatial distribution of the 13 CO 2 discrimination rate over land has considerable impacts on the spatial distribution of the inverted CO 2 sink over land (up to 15 % in some regions), suggesting that reliable models for simulating the spatial distribution of the 13 C discrimination rate over land are needed for effective use of 13 CO 2 data for global carbon cycle inversion.
3. The joint inversion is sensitive to the 13 CO 2 disequilibrium fluxes over both land and ocean.Ignoring these fluxes in the joint inversion causes the inverted total land and ocean sink to increase by 1.18 and 1.96 Pg C year −1 , respectively.The uncertainty in our disequilibrium flux calculation is estimated to be 8.0 and 12.7 Pg C year −1 ‰ for land and ocean, respectively, inducing an uncertainty in the inverted flux of 0.47 Pg C year −1 for land and 0.54 Pg C year −1 for ocean.

Figure 2 .
Figure 2. (a) Gross primary productivity (GPP) distribution in 2003 computed using remote sensing LAI and land cover maps and climate and soil data, and (b) net ecosystem productivity (NEP) distribution in 2003.Both are calculated using the BEPS model.Annual NEP maps from 2000 to 2004 are used to as the prior flux in the inversions.This GPP map is used to distribute the flux uncertainty among the 39 land regions.

Figure 3 .
Figure 3.The annual mean of the total photosynthetic 13 C discrimination ( in Eq. 7) in 2003.

Figure 4 .
Figure 4. Global distribution of the flux-weighted mean age of soil carbon pools (Eq.8).

Figure 6 .
Figure6.Disequilibria between 13 C fluxes to and from the land or ocean surface in 2000.At the land surface, the disequilibrium is the difference between photosynthetic and respiratory discriminations against 13 C, and at the ocean surface it is the difference in 13 C discrimination between the one-way diffusive downward and upward fluxes.

Figure 7 .
Figure 7.Comparison of land and ocean carbon sinks derived from inversions with and without the 13 CO 2 constraint against the Global Carbon Project results (Le Quéré et al., 2013).

Figure 10 .
Figure 10.Comparison between inversion results with and without 13 CO 2 constraint for 21 regions of the globe for the periods of 2002-2004.

Figure 11 .
Figure 11.Comparison between inversion results with and without 13 CO 2 constraint for 30 regions in North America.

Figure 12 .
Figure 12.(a) Difference in ε lph (‰) and (b) the inverted CO 2 flux (g C m −2 year −1 ) between case I and case II, i.e.Case I minus case II.See Sect.2.1.3for the description of these cases.

Figure 14 .
Figure 14.CO 2 concentrations from (a) prior, (b) posterior from the CO 2 -only inversion and (c) posterior from the joint inversion in comparison with observations.The prior concentration is obtained through transport modelling with prior CO 2 fluxes from the terrestrial ecosystems, oceans, fossil-fuel emission and biomass burning.

Figure 15 .
Figure 15.Comparison of prior (a) and posterior (b)13 CO 2 compositions with observations.The prior composition is obtained through transport modelling with prior 13 CO 2 fluxes from the terrestrial ecosystems, oceans, fossil-fuel emission and biomass burning, and the posterior composition is obtained with the CO 2 -13 CO 2 joint inversion (case I).

Table 1 .
Chen et al. (2012)ters are assigned by plant functional types in BEPS.References for the chosen values of these parameters are found inChen et al. (2012).V cmax is the leaf maximum carboxylation rate at 25 • C, J max is the maximum electron transport rate, N is the leaf nitrogen content, χ N is the slope of V cmax variation with N , and m and b are the slope and intercept in the Ball-Berry equation.The peak growing season LAI and clumping index are given as the mean and standard deviation for each plant functional type. *

Table 2 .
Global average ages of soil carbon pools computed by BEPS with consideration of the influences of temperature and soil moisture on the decomposition rates of these pools.

Table 4 .
Comparison of land and ocean disequilibrium coefficients and disequilibrium fluxes calculated in this study with those in previous studies.

Table 5 .
Global isotopic mass budgets averaged for the 2002-2004 period for the prior, double deconvolution, CO 2 -only inversion, and joint inversion (unit: Pg C year −1 ‰).Also shown are ocean and land net fluxes (unit: Pg C year −1 ) for these cases for comparison purposes.For the prior fluxes, the component of each flux are indicated in the parentheses.The isotopic coefficients are same among the cases. .