Interactive comment on “ Regional CO 2 flux estimates for 2009 – 2010 based on GOSAT and ground-based CO 2 observations ” by S . Maksyutov

The authors have conducted an inversion analysis of GOSAT XCO2 and ground-based CO2 data to determine the additional information provided by GOSAT observations, compared to the surface data, to reduce uncertainty in regional CO2 flux estimates. Although GOSAT was launched almost four years ago, issues with biases have made using the data challenging and, as a result, much work is still needed to assess the utility of the data. In this context, the work described in the manuscript is a useful first step in understanding complementarity between GOSAT and existing ground-based CO2 data. I recommend publication of the manuscript in ACP after the authors have adequately addressed the comments below.


Introduction
The recent increase in atmospheric CO 2 concentration is partially abated by carbon uptake by ocean and land, which indicates disequilibrium in CO 2 exchanges between the atmosphere and oceans and between the atmosphere and the terrestrial biosphere (Keeling et al., 1995).The disequilibrium in the terrestrial carbon cycle can be attributed to 20th century warming leading to enhanced nitrogen recycling in the biosphere, coincident with an increase in CO 2 concentrations that facilitates photosynthesis and vegetation functioning (Melillo et al., 2002;Grant et al., 2009).In some regions vegetation recovery is also considered as an important mechanism for net long-term carbon sink (Caspersen et al., 2000;Pacala et al., 2001).The sustainability and amount of the terrestrial sink is still difficult to assess on large regional scales (e.g.Dolman et al., 2012;Gloor et al., 2012) modeling of carbon exchange at the Earth's surface plays an important role in the quantification of the distribution of terrestrial carbon sinks.Analyses of the regional CO 2 fluxes using inverse models of atmospheric transport have proven useful for quantifying the spatial distribution and interannual variability of surface CO 2 fluxes on both global and regional scales (Bousquet et al., 2000;Peters et al., 2007).One particularly important topic under investigation is partitioning the terrestrial carbon sink into (1) sinks in mid and high latitude regions of the Northern Hemisphere, where the warming is most pronounced (Jones and Briffa, 1992), and (2) those in wet tropical and subtropical regions, such as in southern China forests, where net ecosystem production is high (Piao et al., 2009).It was initially found that latitudinal CO 2 gradient suggested stronger Northern Hemisphere sinks (Tans et al., 1990).A more detailed analysis with atmospheric transport and inversion models allocated a large sink to the US (Fan et al., 1998), a result supported by a bottom-up estimate (Pacala et al., 2001).Later, the US sink estimate was reduced in a multi-model inverse modeling study (Gurney et al., 2002) in which a significant net tropical source was also projected.Based on inverse model estimates of the vertical CO 2 gradients, Stephens et al. (2007) supported more moderate estimates of Northern Hemisphere extratropical sinks and near-neutral net CO 2 fluxes in the tropics.However, on a global scale, many gaps appear in remote areas that are not covered by conventional atmospheric CO 2 observation networks (mostly tropical regions), leaving space for large uncertainties in the reconstructed fluxes (Gloor et al., 2000).Rayner and O'Brien (2001) suggested an unconventional solution to the problem of data gaps, arguing that a large number of relatively low-precision satellite observations of atmospheric CO 2 concentration can be used to fill those gaps.This suggestion raised high expectations for the usefulness of remote sensing observations of atmospheric CO 2 .
The launch of the Greenhouse gases Observing SATellite (GOSAT) in 2009, which observes high-resolution spectra of reflected light (Kuze et al., 2009), was followed by continuous efforts to refine retrievals of CO 2 and CH 4 column abundances (Yokota et al., 2009;Bösch et al., 2011;Butz et al., 2011;Yoshida et al., 2011Yoshida et al., , 2013;;O'Dell et al., 2012;Oshchepkov et al., 2012).The availability of GOSAT retrievals validated with Fourier Transform Spectrometer (FTS) observations collected in the Total Carbon Column Observing Network (TCCON) (Morino et al., 2011) provides the scientific community with an opportunity to apply GOSAT data to carbon cycle studies.Theoretical studies on the utility of GOSAT retrievals by Chevallier et al. (2009), Kadygrov et al. (2009), and others suggested that the data can help to fill the gaps in observation coverage if sufficient retrieval accuracy and precision are achieved.
This paper provides an overview of a carbon cycle modeling system that consists of components for modeling atmospheric transport, anthropogenic CO 2 emissions, and terrestrial and oceanic CO 2 exchanges.It further describes the application of the system to the estimation of surface CO 2 fluxes from GOSAT retrievals.The carbon cycle modeling system was developed specifically to utilize the GOSAT data in the inverse model analysis with a delay of a year or less.
In our study, we assessed the utility of the GOSAT X CO 2 retrievals in inverse modeling of surface sources and sinks.We used a recently-improved version of GOSAT X CO 2 retrievals (version 02.00; Yoshida et al., 2013) and the GLOBALVIEW-CO2 ground-based data (2011, hereafter denoted as GV).We estimated monthly regional fluxes and their uncertainties from (1) the ground-based GV data only and (2) both GV and GOSAT X CO 2 retrievals, and compared these two sets of results.The present study used GOSAT data obtained during a one-year period between June 2009 and July 2010, the first year of GOSAT sounding.Section 2 introduces the components of the modeling system.Section 3 briefly describes the GOSAT retrievals and the inverse model.Section 4 presents the results, and Section 5 concludes the paper.

Inverse modeling system components
In this section, we present the components of the inverse modeling system, which are models for simulating (1) atmospheric transport of CO 2 , (2) CO 2 exchange between the atmosphere and oceans, (3) CO 2 exchange between the atmosphere and terrestrial biosphere, and (4) emissions of CO 2 by fossil fuel consumption and cement manufacturing.The a priori flux dataset used in this study is comprised of four components: daily net ecosystem exchange (NEE) predicted by the terrestrial biosphere process model VISIT (Vegetation Integrative SImulator for Trace gases) (Ito, 2010;Saito, M. et al., 2011); monthly ocean-atmosphere CO 2 fluxes generated by an ocean pCO 2 data assimilation system (Valsala and Maksyutov, 2010); monthly CO 2 emissions due to biomass burning stored in the Global Fire Emissions Database (GFED) version 3.1 (van der Werf et al., 2010); and monthly fossil fuel CO 2 emissions obtained via combining the high-resolution Open source Data Inventory of Anthropogenic CO 2 emission (ODIAC) dataset (Oda and Maksyutov, 2011) and the Carbon Dioxide Information Analysis Center's (CDIAC) monthly 1 • × 1 • resolution dataset (Andres et al., 1996(Andres et al., , 2011)).Each of these component flux datasets was prepared specifically for this 2009-2010 analysis period.

Model of the carbon cycling in the terrestrial biosphere
VISIT is a prognostic biosphere model (Ito, 2010;Saito, M. et al., 2011) that simulates carbon exchanges between the atmosphere and biosphere and among the carbon pools within terrestrial ecosystems at a daily time step.The carbon pools in the model consist of five compartments: foliage, stem and branch, root, litter, and soil.Modeling of plant CO 2 assimilation in VISIT is based on a model of light extinction in the canopy, following the formulation of Monsi and Saeki (1953).Maximum photosynthetic uptake rate is influenced by temperature, atmospheric CO 2 concentration, and soil moisture.Autotrophic respiration is formulated as the sum of growth respiration and maintenance respiration.Growth respiration is simulated as the cost to produce new biomass, while maintenance respiration is represented as a function of ground surface temperature.Heterotrophic respiration is the sum of respiration from two layers, litter and humus, which is regulated by soil temperature and soil moisture at each depth.Litterfall from foliage, stems and branches, and roots is calculated by a simple parameterization on the basis of the carbon mass of each component.NEE, which is one of the a priori fluxes required in forward CO 2 concentration simulations with an atmospheric transport model, is given as the difference between ecosystem respiration and gross primary productivity.Here, ecosystem respiration is the sum of autotrophic respiration and heterotrophic respiration.A positive value of NEE indicates CO 2 release to the atmosphere from the terrestrial biosphere, whereas a negative value indicates CO 2 uptake from the atmosphere.
VISIT is driven by reanalysis/assimilation climate datasets provided by the Japan meteorological Agency (JMA): the Japan 25-year reanalysis (JRA-25)/JMA Climate Data Assimilation System (JCDAS) (Onogi et al., 2007) for the period 1979 -present.The meteorological data that drives VISIT include downward shortwave radiation at the surface, total cloudiness, 2 m air temperature, ground surface temperature, soil temperature at depths of 10 cm and 200 cm, specific humidity, precipitation, and wind velocity.The JRA-25/JCDAS data are provided at a T106 spatial resolution at 6 h temporal resolution.All of the JRA-25/JCDAS data were converted to daily mean values at a 0.5 • × 0.5 • grid resolution using an interpolation, and then used as forcing data for VISIT.Biases in JRA-25/JCDAS precipitation data were corrected following the method of Saito, M. et al. (2011).The model was initially run for a spin-up of approximately 2000 yr to reach equilibrium in the carbon pools, by repeating JRA-25/JCDAS forcing with varying atmospheric CO 2 .Then the daily physiological processes were simulated for the period starting in 1979.Global vegetation was classified into 16 plant functional types at a 0.5 • × 0.5 • grid resolution using the Moderate Resolution Imaging Spectroradiometer (MODIS) land cover product (Friedl et al., 2002).At each model grid, all physiological processes involve the effect of vegetation fractional coverage up to the fourth dominant biome.
We used an optimized VISIT model to simulate daily NEE variability.The optimization method and results are described in detail by Saito et al. (2013), and are presented only briefly here.Physiological parameters of the model were optimized to fit the observational data using a Bayesian inversion approach, which is an extension of a method by Nakatsuka and Maksyutov (2009).Thirteen selected physiological parameters were optimized for 16 plant functional types.The application of GV CO 2 seasonal cycle as the only constraint led to low biased net primary productivity (NPP).Accordingly, biomass and NPP observations had to be included as additional constraints in order to reproduce the seasonal cycle of CO 2 while keeping NPP and biomass within observed parameter ranges.Data of seasonal variability in atmospheric CO 2 concentrations, annual mean aboveground biomass (AGB), and NPP were assimilated into VISIT.These data were derived from GV, International Institute for Applied Systems Analysis (IIASA) global biomass map (Kindermann et al., 2008), and Global Primary Production Data Initiative (GPPDI) (Scurlock et al., 1999;Olson et al., 2001), respectively.An atmospheric tracer transport model (Maksyutov et al., 2008) was used for the computation of atmospheric CO 2 variability.Both VISIT and the atmospheric transport model were run on 2.5 • × 2.5 • grid resolution for the optimization.The misfits between the simulated and observed monthly mean CO 2 concentrations and other parameters (NPP and AGB) were minimized iteratively.In each iteration, the Jacobian matrix was first estimated by calculating the sensitivity of the simulated monthly concentrations, NPP and AGB to small changes in the VISIT parameters.Then an optimal set of parameters was found by solving a linear inverse problem.The iterations were repeated several times because the VISIT parameters and the modelsimulated CO 2 fluxes are not linearly related.
The VISIT model is under continuing development and validation against observations.Biomass amounts, plant productivities and the temporal variations of NEE simulated by VISIT were discussed in multi-year and multimodel comparison studies by Ito and Sasai (2006), Ito et al. (2010), Ichii et al. (2010), Ichii et al. (2013), Piao et al. (2012).Saeki et al. (2013), in their forward transport and inverse modeling study, demonstrated that the seasonal cycles of CO 2 concentration simulated with VISIT NEE at selected mid to high latitude sites and those of observed values were in reasonable agreement.The modeled amplitudes of the seasonal cycles at some high latitude sites were, however, shown to be smaller than the observed.Valsala et al. (2013) compared the intraseasonal variability of NEE simulated by VISIT and Carbon-Tracker (Peters et al., 2007) over Indian monsoon season and found that the two estimates were similar in intra-seasonal variability but different in seasonality.
For the analysis period from 2009/06 to 2010/05 the annual net flux was −0.7 GtC yr −1 .The time series of model simulated CO 2 concentration at Mauna Loa Observatory, based on VISIT NEE, is compared with corresponding GLOBALVIEW data in Fig. 1.Also shown in the figure are biomass maps simulated by VISIT and estimated by IIASA.

Variational assimilation system for simulating the global pCO 2 maps and surface ocean-atmosphere fluxes of carbon
The magnitude of annual atmospheric CO 2 flux into the oceans is estimated to be 1.5 to 2 PgC yr −1 (Gurney et al., 2004, Gruber et al., 2009).Therefore, considerable efforts have been given to the preparation of the oceanic fluxes used in this study.
The air-sea CO 2 flux component used in this study was taken from an optimal estimate of oceanic CO 2 fluxes derived from the original work of Valsala and Maksyutov (2010).This dataset was produced by simulating dissolved inorganic carbon (DIC) with a simple ocean biogeochemical model and constraining DIC to observations through a variational assimilation method.The data is available from 1996 to near real-time.The essential components of the oceanic model are described in below.
In the work of Valsala and Maksyutov (2010), a simple offline ocean tracer transport model (OTTM) described by Valsala et al. (2008) was coupled with a simple onecomponent ecosystem model based on phosphate cycling (McKinley et al., 2004) and an abiotic carbon cycle model of OCMIP-II (Ocean Carbon Cycle Intercomparison Project, Orr et al., 1999), and was used to simulate the air-to-sea CO 2 fluxes.The model surface ocean DIC values were then constrained with corresponding observational values that were derived from the observed partial pressure of surface ocean CO 2 (pCO 2 ) obtained via numerous ship-underway sampling summarized in Takahashi et al. (2011) database.The assimilation scheme, which consists of a variational method derived from the work of Ikeda and Sasai (2002), minimizes model-observation differences in surface ocean DIC (or pCO 2 ).The transport model was run with offline ocean currents provided by GODAS re-analysis data products (Behringer and Xue, 2004).The offline data fed into the system are ocean current velocities, temperature, salinity and other physical parameters that were derived from the re-analysis data at a five-day time interval.The OTTM tracer transport model was run at 1 • × 1 • resolution and 40 vertical levels.The first 26 levels are in upper 300 m of the ocean.The use of the offline re-analysis input fields for running the transport model enabled us to simulate the air-sea CO 2 flux in near real-time.The simulated ocean DIC values were then corrected with observational data of surface ocean pCO 2 using a two-way constraining process in the assimilation.The model surface ocean pCO 2 are constrained strongly whenever the ship-track underway sampling is available.In addition to this rather "strong" constraint, the climatological maps of monthly mean pCO 2 derived from Takahashi et al. (2009) were also used to constrain the surface ocean pCO 2 as a "weak" constraint.This two-way correction applied to model surface ocean pCO 2 (i.e.effectively to the DIC) reduces the model biases as well as DIC errors.60% of the annual mean model biases were eliminated in the assimilation, and 40-60 % of the cumulative seasonal errors were also reduced at regional scales (see also Valsala and Maksyutov, 2010).Valsala and Maksyutov (2010) used offline data from the Geophysical Fluid Dynamics Laboratory (GFDL) re-analysis products to drive the transport model.However, because the GODAS re-analysis dataset is updated and made available in near real-time, we chose to use this data set for our operational estimates of air-sea CO 2 fluxes.Two main advantages of employing the OTTM-derived optimal air-sea CO 2 fluxes can be summarized as follows: (1) OTTM-derived fluxes are monthly varying and available in near real-time, and (2) the air-sea CO 2 fluxes thus have signatures of interannual variability, as compared to the monthly climatological maps of air-sea CO 2 fluxes based on Takahashi et  2009), which are often employed in inversion setups (e.g.Gurney et al., 2004).
The interannual variability of the ocean-atmosphere exchange simulated with OTTM has been analyzed by Valsala et al. (2012a), who found a persistent quadra-pole pattern of CO 2 flux interannual variability (IAV) in the North Pacific varying at Pacific Decadal Oscillation scale.Valsala and Maksyutov (2013) and Valsala et al. (2012b) analyzed the flux IAV in the northern Indian ocean.Ishii et al. (2013) compared the assimilated flux dataset with multiple bottom-up and inverse modeling estimates within a framework of REC-CAP (REgional Carbon Cycle Assesment and Processes) project (Canadell et al., 2011), and showed that the interannual variation of the tropical Pacific fluxes correlates with atmospheric inverse model estimates.
Figure 2 shows the global average of air-to-sea CO 2 fluxes from June 2009 to May 2010, global integrated CO 2 sink and data uncertainties for individual months, used in this study.An annual mean of 2.02 PgC yr −1 of CO 2 sink is resolved in the optimized flux for the period of inversion.

Emissions dataset for fossil fuel CO 2 emissions
Similar to many other inversion studies (e.g.Gurney et al. 2002), fossil fuel CO 2 emissions (emissions due to combustion of fossil fuels and cement manufacturing) in our inverse estimation are not solved, but rather prescribed.Thus, fossil fuel CO 2 emissions need to be accurately given for better flux estimation (Gurney et al. 2005).
To prescribe fossil fuel CO 2 emissions, we used an updated version of the ODIAC dataset (Open-source Data Inventory for Anthropogenic CO 2 ; Oda and Maksyutov ( 2011)) prepared at 1 × 1 degree resolution on a monthly basis.Monthly estimates of national total emissions are available from the Carbon Dioxide Information and Analysis Center (CDIAC) of the U.S. Department of Energy (Boden et al., 2011; http://cdiac.ornl.gov/trends/emis/overview2008.html, last access: Aug 2, 2012).These emissions estimates are projected up to year 2010 using British Petroleum statistical data (British Petroleum p.l.c., 2011) and CDIAC's preliminary estimate (http://cdiac.ornl.gov/trends/emis/prelim2009 2010 estimates.html,last access: 2 August 2012).The CDIAC estimates comprise of emissions from several categories: fuels (solid, liquid, and gas), cement production, gas flaring, and international bunkers.Emissions from solid, liquid, and gas fuels and cement production were then spatially distributed using power plants profiles (geographical location and emissions) given by CARMA database (CARbon Monitoring and Action, www.carma.org,last access: 6 August 2012) and satellite-observed nightlight data collected by U. S. Air force Defense Meteorological Satellite Project (DMSP) satellites (Elvidge et al. 1999).The nightlight data were processed by National Oceanic and Atmosphere Administration (NOAA) National Geophysical Data Center (NGDC) (Ziskin et al, 2010).This distribution method, compared to previous studies such as Andres et al. (1996), allows us to map emissions over land at a high spatial resolution (up to 1km) (Oda and Maksyutov, 2011).At spatial resolutions of global transport simulation, the resulting spatial distribution agrees well with that of North American emissions data product Vulcan (Gurney et al., 2009) (Oda and Maksyutov, 2011).
In addition to spatial distribution, seasonality in land emissions (solid, liquid, gas and cement) was adopted from CDIAC's monthly 1 • × 1 • emissions dataset (Andres et al., 2011; data available at: http://cdiac.ornl.gov/epubs/fossilfuel CO2 emissions gridded monthly v2011.htmllast access: 2 August 2012).Similarly, AERO2k inventory was used for emissions from aviation.As a simplification, aviation emissions are treated as surface emissions in the transport model, as it is a minor component of the total emissions.Gas flaring and marine bunker emissions do not have seasonality in ODIAC dataset.
A global picture of fossil fuel CO 2 emissions is shown in Fig. 3.The figure was drawn using emission data reduced to 5 km × 5 km (2.5 arc min) resolution.Similarly to the emission map by (Oda and Maksyutov, 2011) intensive emissions are found in northern hemisphere especially over industrialized countries and regions.

Emissions of CO 2 by biomass burning and forest fires
Large contribution of fire processes to the interannual variability of the global carbon cycle is known from high correlation between inversion-reconstructed flux anomalies and satellite based fire estimates (Patra et al., 2005a).Satellite based estimates of carbon emissions due to forest fire and biomass burning are provided by Global Fire Emissions Database (GFED 3.1) as described in van der Werf et al. ( 2010) and Giglio et al. (2010).GFED 3.1 used a combination of active fire observations from multiple satellites, 500 m MODIS burned area maps, local regression, and regional regression trees to produce a hybrid, global, monthly burned area data set from July 1996 to December 2010.Annual totals derived from these data show good agreement with independent annual estimates available for Canada and the United States (Giglio et al., 2010), and Russia (Shvidenko et al, 2011).The global annual burned area for the period 1997-2008 varied between 330 and 431 Mha, with a maximum occurring in 1998 and the minimum in 2008.The most extensive burning consistently occurred in Africa, with an average of 250 Mha burned on the continent each year.This represents about 70 % of the total area burned worldwide annually.The lowest interannual variability in the extent of burned areas occurred in the African savannas.Regions of high interannual variability included Australia, the United States, and the boreal forests of both Asia and North America.Burned area maps were produced from the 500-m MODIS atmospherically-corrected Level 2G surface reflectance product (Vermote et al., 2002), the MODIS Level 3 daily active fire products (Justice et al., 2002) and the MODIS Level 3 96-day land cover prod-uct (Friedl et al., 2002), by finally applying the Giglio et al. (2009) burned area mapping algorithm.Algorithm validation for 500-m burned area maps is completed for Southern Africa, Siberia, and the Western United States through comparison with Landsat imagery (Giglio et al., 2009).For active fire counts a MODIS monthly Climate Modeling Grid (CMG) fire product at 0.5 deg spatial resolution was used.Detection of fire and burnt area, as well as fraction of the biomass and biomass debris is accompanied by relatively large uncertainties (Giglio et al., 2010) and contribute to the overall uncertainty of the inverse model estimates.

Atmospheric tracer transport model
We used the National Institute for Environmental Studies global atmospheric tracer Transport Model (NIES-TM) to run forward simulations of atmospheric CO 2 for the inverse modeling of surface CO 2 fluxes.NIES-TM is an off-line model driven by reanalysis data, which combines JRA-25 for the period 1979-2004 and a real-time operational analysis by the JMA Climate Data Assimilation System (JC-DAS) for the period 2005-present (Onogi et al., 2007).The JRA-25/JCDAS data used in our model is provided on T106 Gaussian horizontal grid (320 × 160 grid points) with 40 hybrid σ -p levels and the 6-hour time step.In the version used in this study (version NIES-08.1i), a flexible hybrid sigma-isentropic (σ -θ) vertical coordinate system was implemented, which combines terrain following vertical coordinate in the troposphere and isentropic vertical coordinate in the stratosphere above the level of 350K to ensure that isentropic surfaces do not intersect the Earth's surface (Belikov et al., 2013a).
A scheme based on radiative heating and cooling was implemented to calculate vertical transport in the stratosphere, because such a scheme produces a better representation of the meridional circulation for long-term simulations, as compared with use of vertical winds from reanalysis (Weaver et al., 1993).Air-ascending rates controlled by climatological heating rate derived from JRA-25/JCDAS reanalysis data were adjusted to fit to the observed mean age of the air in the stratosphere.The model employs a reduced latitudelongitude grid scheme in which the grid size is doubled approaching the poles.This approach overcomes the Courant condition limitation on a model time step in the polar regions, caused by small grid size associated with meridional convergence, maintains a reasonably high integration timestep and ensures adequate model performance (Belikov et al., 2011).The model uses a flux-form advection algorithm with a second-order van Leer scheme.To obtain mass conservation in the numerical scheme, the horizontal mass fluxes derived from the meteorological dataset are balanced with the surface pressure tendency using the method developed by Heimann and Keeling (1989).As in previous model versions (Maksyutov et al., 2008;Belikov et al., 2011) the parameterization of the turbulent diffusivity follows the approach used by Hack et al. (1993).Under the assumption that the planetary boundary layer (PBL) is well mixed the turbulent diffusivity is set to a constant value of 40 m 2 s −1 in the PBL.The free-tropospheric diffusivity is estimated as a function of the Richardson number.To separate the transport processes in the PBL and above it we used 3-hourly PBL height data obtained from the European Centre for Medium-Range Weather Forecasts (ECMWF) Interim Reanalysis dataset (Dee et al., 2011).
We implemented the Kuo-type penetrative cloud convection scheme proposed by Tiedtke (1989), with entrainment and detrainment processes on convective updrafts and downdrafts to simulate deep convection (Belikov et al., 2013b).Calculation of cumulus mass-flux is based on the method developed by Austin and Houze (1973), in which the amount of air lifting in an updraft core of a cumulus cell is related to the amount of precipitation the cumulus cell produces.The mass of air transported upward within the cells was computed from the conservation of moisture, using the convective precipitation rate from the JRA-25/JCDAS dataset.
The model was evaluated against GLOBALVIEW-CO2, GLOBALVIEW-CH4, World Data Centre for Greenhouse Gases (WDCGG), balloon-borne and aircraft observation data (Belikov et al., 2011(Belikov et al., , 2013a)), as well as through the Comprehensive Observation Network for Trace gases by AIrLiner (CONTRAIL) (Niwa et al., 2011) and TransCom-CH 4 transport model intercomparison (TMI) studies (Patra et al., 2011, Belikov et al., 2013b).Implementation of hybrid sigma-isentropic vertical coordinates with a radiative balance scheme for vertical transport allows simulation of vertical profiles and vertical propagation of seasonal variations of tracers in the free troposphere and in the lower stratosphere which appears to be in good agreement with aircraft and balloon-borne observations.In general NIES-TM performance is consistent with the TransCom-CH 4 and CON-TRAIL intercomparison participating models (Niwa et al., 2011;Patra et al., 2011).Comparisons with balloon-borne observations over Sanriku, Japan in 2000-2007 revealed that the tracer transport simulations in the upper troposphere and lower stratosphere are performed with accuracies of ∼ 5 % for CH 4 and SF 6 , and ∼ 1 % for CO 2 compared with the observed volume-mixing ratios (Belikov et al., 2013a).
The model is able to reproduce the seasonal variations in CO 2 and CH 4 surface concentrations.More accurate results were obtained for CH 4 at regions located some distance away from multiple emission sources.For other sites, where high emissions and local meteorology play a major role, it proved difficult to reproduce the CH 4 surface concentrations, especially in winter, which indicates excessive near-surface vertical mixing under stable conditions.
Modeled dry-air column-averaged CO 2 and CH 4 values (X CO 2 and X CH 4 ) were also evaluated against daily groundbased high-resolution FTS observations measured at 12 sites of the Total Carbon Column Observing Network (TCCON) for the period from January 2009 to January 2011.Modeled data convolved with scene-dependent averaging kernels were able to reproduce the seasonal and inter-annual variability of X CO 2 and X CH 4 with correlation coefficients of 0.8-0.9 and 0.4-0.8, and biases ±0.2 % and ±0.5 % (without Sodankyla's data), respectively (Belikov et al., 2013a).
Modelled X CO 2 and X CH 4 were also evaluated against daily ground-based high-resolution FTS observations measured at 12 sites of TCCON for the period from January 2009 to January 2011.Modeled data convolved with scenedependent averaging kernels were able to reproduce the seasonal and inter-annual variability of X CO 2 and X CH 4 with correlation coefficients of 0.8-0.9 and 0.4-0.8, and biases ±0.2 % and ±0.5 % (without Sodankyla's data), respectively (Belikov et al., 2013a).
The atmospheric transport was simulated at resolution of 2.5 • × 2.5 • on 32 vertical levels.The model can use surface fluxes at hourly, daily and monthly temporal resolutions and 1 • × 1 • spatially.The fluxes simulated by surface flux models were converted to 1 • × 1 • resolution when necessary.The OTTM oceanic fluxes come as monthly mean data.VISIT model was run at 0.5 • × 0.5 • resolution at a daily time step.The VISIT fluxes were converted to 1 • × 1 • daily fields.Use of daily mean fluxes as a substitute for diurnally varying fluxes is supported by estimates by Olsen and Randerson (2004), who found that early afternoon satellite observations of X CO 2 represent well the daily mean fluxes.GFED and ODIAC fluxes were provided to the transport model at 1 • × 1 • resolution and a monthly time step.The 1 • × 1 • fluxes were remapped to 2.5 • × 2.5 • fields inside the transport model.

Inverse modeling scheme
The theoretical basis for the technique is Bayes' theorem (see e.g.Tarantola, 2005), with which the "optimal" or a posteriori state of a set of parameters is deduced from a priori knowledge about those parameters and measured data values.For the case of estimating surface fluxes of CO 2 , which is approximated to be chemically inert, the relationship between the measured data values and their theoretical predictions based on physical process modeling is linear.The relationship can be expressed in matrix form as where d obs is the concentration vector recorded at measurement locations, m denotes modeled source strengths in predefined regions, and is a matrix that maps the source strength field onto that of concentrations.The elements of matrix G are given as changes in concentrations at each of measurement sites with respect to unit pulse emissions from each of the pre-defined regions.These elements are obtained by running forward a set of unit pulse emissions with an atmospheric tracer transport model (e.g.Rayner et al, 1999).
Vector v is the misfit between the observations and the model-predictions, which is comprised of measurement uncertainty and error in the simulation of atmospheric transport.The aim here is to find m that best describes d obs .In the context of Bayesian inference, where Gaussian probabilities are assumed for measurements and a priori parameters, the measure of the fit between modeled source strengths m and measurements is expressed as a cost function L(m) (also see Tarantola, 2005): where m prior denotes the vector of the a priori source strengths, and C D and C M are the error covariance matrices of the measurements and the a priori source strengths, respectively.The optimal state of the modeled source strengths, m , exists at the minimum of this measure.Taking the derivative of L with respect to m and setting it to zero yields m = m prior (3) Further, differentiating L with respect to m for the second time gives the a posteriori error covariance, C M , which is a measure of the convexity of the function L.
The values assigned to the elements of vector d obs in the present study were monthly-mean, surface-based GV data and bias-corrected GOSAT X CO 2 retrievals gridded to 5 • × 5 • cells and averaged on a monthly basis.The GV data are a product generated with a technique developed by Masarie and Tans (1995), which incorporates interpolated and/or extrapolated values with measurements such that the resulting smoothed concentration time series become seamless in time.The reason behind the choice of monthly mean GV data, instead of using simple averages of available flask observations in each month, as in a study by Rödenbeck et al. (2003), is to minimize the impact of temporal discontinuities among those observations on the flux estimation.The model simulated concentration at each observation location of GV and GOSAT X CO 2 was obtained by performing linear interpolation, in space and time, of the model-predicted concentration field (updated at a time step of 10-15 min).Cells with three or more X CO 2 retrievals per month were selected.Prior to monthly averaging, large GOSAT X CO 2 outliers were removed via comparisons with climatological X CO 2 values, derived from an ensemble of forward simulation results by six different transport models (Gap-filled and Ensemble Climatology Mean: GECM) that was nudged to surface-based observations (Saito, R. et al., 2011).The observation errors for the monthly mean X CO 2 retrievals, specified in the diagonal elements of matrix C D , were determined as the standard deviations of GOSAT X CO 2 retrievals found in each of the 5 • × 5 • grid cells in a month.We took account of errors associated with the retrieval of X CO 2 values and the forward atmospheric transport simulation by setting the minimum of the observation error for GOSAT X CO 2 retrievals at 3 ppm, which consists of an uncertainty associated with the retrieval of GOSAT X CO 2 (2 ppm) and that of forward X CO 2 modeling (1 ppm) The GV sites were selected by comparing GV data against concentrations predicted by NIES-TM over the analysis period.We picked the sites whose RMS modelobservation misfits were less than 2 ppm.As an observation error estimate, the GV residual standard deviation (stored in the GV dataset) was assigned to each of the selected sites.We gave less weight at GV sites whose observational record completeness was less than 70 % by tripling their data errors.The minimum error for the GV data was set at 0.3 ppm.Altogether, 220 GV data time series were selected for this estimation.Some sites such as aircraft observation sites contain several time series per site.
The diagonal elements of the matrix C M were prescribed as follows.The uncertainty of the terrestrial a priori flux was set at twice the standard deviation of the VISIT model monthly NEE (1 • × 1 • resolution) values for the past 30 yr.The uncertainty of the oceanic a priori flux was determined as the standard deviation of the OTTM-assimilated oceanic flux (1 • × 1 • resolution) for the period 2001-2009, and the climatological mean data of Takahashi et al. (2009).
In the TransCom 3 CO 2 inversion intercomparison, Gurney et al. ( 2003) assigned growing season net fluxes (GSNF; the sum of monthly-mean exchanges for months exhibiting net uptake) as terrestrial prior flux uncertainties (GSNF were based on NEE predicted by the CASA model).The reason behind it was that GSNF provide ecologically relevant upper bounds for annual-mean terrestrial flux.For oceanic fluxes, Gurney et al. (2003) set the uncertainties at 140 % of the climatological net oceanic exchanges, which are approximately double the amount suggested by Takahashi et al. (2002).Our approach of using standard deviations of VISIT NEE and OTTM oceanic fluxes is similar to their case in finding reasonable upper limits of naturally varying fluxes and assigning them as boundaries in the flux estimation.Our boundaries reflect natural variability in the past several decades (30 yr for terrestrial biosphere and 10 yr for ocean).
The size of vector m in the present study was set to the number of source regions (64 regions) times the number of analysis months (14 months).The 64 source regions consist of 42 subcontinental-scale terrestrial regions and 22 ocean basins (Patra et al. 2005b).The boundaries of these source regions are shown in Fig. 4. The dimension of matrix G is then determined as the size of vector m multiplied by that of vector d obs .For implementing matrix operations involved in Eq. ( 3) efficiently, we employed a variant of the fixed-lag Kalman Smoother scheme (FLKS) formulated by Bruhwiler et al. (2005).The basis for this scheme is the fact that in atmospheric tracer transport simulations, the signals of unit pulse emissions detected at measurement sites decay rapidly within the first few months and are blended into the background state thereafter.The idea is to obtain a posteriori fluxes via estimating m incrementally with a subset of G and d obs in a specified time-window.Using the FLKS setup with the same 64 source region boundaries, Koyama et al. (2009) evaluated the influence that differences in the length of the time window have on a posteriori monthly flux estimates.Comparing results obtained using window lengths of 1 to 6 months, they concluded that a posteriori fluxes and their uncertainties estimated with three-month or longer windows converged quite strongly; Bruhwiler et al. (2005) arrived at a similar conclusion.Based on these findings, we chose to use a three-month duration of monthly observations and atmospheric transport simulations in each time window.
In implementing the scheme, only random errors associated with the prior fluxes and observations, as specified in matrices C D and C M , are accounted for.Thus, biases in the observations need to be removed prior to the optimization.We placed two additional columns in matrix G that correspond to global CO 2 offsets of (1) the initial model concentrations with respect to GV data, and (2) the initial model concentrations with respect to GOSAT X CO 2 retrievals.In an initial optimization run, these two offsets were determined separately every month.The averages of the monthly offset values were then used in the subsequent run to remove the global offsets.

GOSAT X CO 2 retrievals
The main observational instrument aboard GOSAT is the Earthward-looking TANSO-FTS that measures surfacereflected sunlight and emitted thermal infrared radiation at wavelengths in the range 0.76-14.3µm.The design and functions of the instrument are described in detail by Kuze et al. (2009).Sampled spectra recorded in the 0.76 µm oxygen absorption band and the 1.61 µm CO 2 absorption band were used in an earlier version of the NIES Level 2 operational retrieval algorithm (version 01; described by Yoshida et al., 2011) to retrieve X CO 2 global distributions.The retrieved X CO 2 values exhibited promising characteristics, including distinct north-south gradients and seasonal variability.However, the X CO 2 values were found to contain a significant negative bias of 8.85 ±4.75 ppm (Morino et al., 2011) as compared with reference data collected at TCCON sites (Wunch et al., 2011a), where sun-viewing high-resolution FTSs are installed.Later, Uchino et al. (2012), using their lidar observations of aerosol particles, showed that assumptions made in version 01 of the retrieval algorithm on the vertical distributions of thin cirrus and aerosols are oversimplified, thereby contributing to the large bias.They proved that the issue could be mitigated significantly by the use of aerosol/cirrus optical properties retrieved simultaneously with spectra in the 2.06 µm band.Further, through investigating GOSAT spectra sampled over 2.5 yr, Yoshida et al. (2012) discovered a time-dependent degradation of TANSO FTS's radiometric accuracy, which they successfully modeled for use in the retrieval algorithm implementation.These new findings, along with other improvements, were incorporated into the NIES Level 2 operational retrieval algorithm.The updated Level 2 X CO 2 retrievals (version 02.00), processed from an improved GOSAT spectral dataset (Level 1B data, version 141.141, covering 14 months from June 2009 to July 2010) were shown to have a much smaller bias of −1.20 ± 1.97 ppm.However, the causes of the remaining bias require further investigation.Wunch et al. (2011b) made an attempt to assess and correct spatially-and temporally-varying biases in GOSAT X CO 2 retrievals using an empirical regression model with which they correlated spurious variabilities in X CO 2 retrievals with surface albedo, difference between the analyzed and retrieved surface pressure, airmass, and the oxygen A-band spectral radiance.A similar analysis is being performed on the GOSAT Level 2 X CO 2 retrievals used in this study.The corrections obtained through the regression analysis, however, still need improvements before being applied to inverse modeling since the derived parameter regression slopes have large uncertainty due to the fact that the number of GOSAT retrievals that match TCCON data is limited.From the view point of global-scale flux estimation, it is critical and highly desired that biases coming from the dependence of X CO 2 retrievals on airmass characteristics, which can change with time and location over the globe, are well-compensated.More detailed regression analyses may be possible by using either data from the newly established TCCON sites and/or well characterized model-simulated CO 2 fields.In addition to the ongoing bias analysis, efforts are also being made at reducing the bias itself through improving the X CO 2 retrieval algorithms.The outcome of these efforts will be reflected in the future updates of the X CO 2 retrieval dataset.For this study, we therefore corrected the bias by raising each X CO 2 value by the global mean GOSAT-TCCON difference of 1.20 ppm prior to the use in inverse modeling, assuming that the bias is uniform throughout the globe and the observation period.
Figure 5 shows the number of GOSAT X CO 2 retrievals per each of 5 The distribution of the data number density changes with season owing to the occurrence of clear sky days and local solar zenith angle that determines the northern-and southern-most bounds of the GOSAT measurement.Note here that regions above ∼ 50 • N latitude (the northern parts of North America and Eurasia) during fall and winter months saw very small numbers of GOSAT retrievals therefore the flux inference for those regions during these months was reliant on the GV data.Figure 6 displays GOSAT X CO 2 retrievals in the form of input to the inverse modeling scheme (gridded to 5 • × 5 • cells and averaged on a monthly time scale).Only the cells with three or more retrievals per month are shown here.The monthly mean GV values are also shown in the figure in circles.Prior to monthly averaging, we corrected the bias in the GOSAT X CO 2 retrievals by raising each X CO 2 value by 1.20 ppm, assuming that the bias is uniform throughout the globe and the observation period.Evaluating the reasonableness of this assumption is a subject of ongoing studies.

Treatment of GOSAT averaging kernel
To account for the vertical sensitivity of the GOSAT measurement in the inverse modeling, we applied the averaging kernel, derived in the retrieval of X CO 2 , to each of the vertical concentration profiles simulated with NIES-TM.As was described by Connor et al. (2008), a model-simulated X CO 2 concentration X m CO 2 , which reflects the measurement vertical sensitivity, is given as: where X a CO 2 denotes a priori X CO 2 values defined in the X CO 2 retrieval, A is a matrix containing the CO 2 elements of the averaging kernel, x m and x a denote the elements of the modeled and the a priori vertical CO 2 profile, respectively.h is the pressure weighting function, a vector containing the dry air partial column abundance of each retrieval layer normalized to the total dry air column abundance.The calculation of the pressure weighting function is described in Appendix B of a report by Yoshida et al. (2009).

Results and discussion
Using the 14-month-long GOSAT Level 2 X CO 2 retrievals (version 02.00) and the GV data in the 3-month-window FLKS scheme, we inferred monthly fluxes for the 64 subcontinental regions for 12 months between June 2009 and May 2010.The forward concentration simulation was initialized with a GECM-derived 2-D global concentration field, and then pre-run was performed for 3 months before the start of the inversion period.A total of 9106 observations were available for estimating 768 monthly fluxes (64 regions ×14 months), of which 6125 were gridded monthly-mean GOSAT X CO 2 retrievals and 2981 were monthly-mean GV data.
The reduction in the a priori flux uncertainty corresponds to the degree to which observations used in the inference contributed to constraining the surface fluxes.The reduction is often expressed by contrasting the diagonal parts of the a posteriori error covariance matrix, C M , to that of the a priori one, C M .Here, we rather chose to consider the uncertainty reduction attained by the addition of the GOSAT X CO 2 retrievals to the GV data.Following Rayner and O'Brian (2001) and Takagi et al. (2011), we express the uncertainty reduction (UR) as: where the units of UR are in %, and σ GV and σ GV + GOSAT denote the uncertainties in the monthly fluxes estimated from the GV data only and those from both the GV data and the GOSAT retrievals, respectively.For this evaluation, we implemented the inversion scheme using only the GV data to obtain flux estimates and the value of σ GV .Figure 7  (3) the size of C D , which reflects the number of observations available for constraining the fluxes.Note that in the current inversion setup the uncertainties specified for GV data and that for GOSAT retrievals can differ by as much as one order of magnitude (e.g. the minimum uncertainty set for GV data and GOSAT retrievals is 0.3 and 3.0 ppm, respectively).This implies that the GV data have much greater weight in constraining regional fluxes.Also, we consider that there is approximately one-order-of-magnitude difference between the uncertainties prescribed to land and ocean fluxes.These differences contribute to creating strong region-to-region or land-to-ocean contrasts in UR values, as observed in Fig. 7. Regions that are far from ground-based observation networks but are covered by GOSAT retrievals (e.g.regions 29 and 17; see Fig. 4 for identifying the regions) show higher UR values, with a maximum UR of 61 % for region 29 in October 2009 (not shown in the figure).However, the UR values for the North American and Australian regions (regions 5-8 and 35-38) barely exceed ∼ 15 %, despite the fact that GOSAT retrievals were constantly available within and around these regions throughout the 1-year analysis period (see Fig. 6).This represents a case in which the constraint provided by the GV data prevails over that provided by the GOSAT X CO 2 retrievals.Thus, higher URs in the figure highlight regions whose a posteriori fluxes were constrained by the GOSAT retrievals more strictly than those in other regions (Middle East, Asia, Africa, and South America).In light of the GOSAT mission objectives, Fig. 7 indicates what the satellite was designed to perform in complementing the groundbased observations.However, care must be taken in evaluating the flux values, as these remote regions coincide with locations where the validation of GOSAT retrievals is not currently possible and the retrieval of X CO 2 values itself is challenged by higher local surface albedo and/or contamination by aerosols.
Figure 8 shows the monthly a posteriori fluxes described above.The quantity presented here is the sum of a priori fluxes (terrestrial biosphere exchange or ocean exchange + anthropogenic emissions + forest fire emissions) and the correction to the a priori flux determined via the optimization.The net influence of the addition of the GOSAT X CO 2 retrievals to the GV data on the flux estimation is made visible by taking the difference between the flux correction obtained with GV data only and that with both the GV data and the GOSAT retrievals (see Fig. 9).As indicated by the distribution of UR rates shown in Fig. 7, the major changes in the correction fluxes occurred in the data-poor regions of surface observation networks.Among the changes, those associated with very low URs (< ∼ 15 %; e.g.tropical America, central Africa, and southern Asia in summer and spring) need to be isolated as they are most likely the result of flux balancing or compensation that took place through the optimization of neighboring regional fluxes.The remaining changes, which are associated with higher URs (e.g.northern Eurasia during summer, the Middle East, southern Africa, and central Asia) and therefore more strongly linked to the GOSAT X CO 2 retrievals, are explored further here.The site-by-site GOSAT data validation activities showed that the X CO 2 retrievals now agree reasonably well with the TCCON reference dataset, with a global mean bias of −1.20 ± 1.97 ppm.Evaluating the validity of spatiotemporal changes in the global distribution of GOSAT-based X CO 2 concentrations will, however, continue to present challenges, perhaps until the cross-validation of similar space-based CO 2 measurements is possible.Here, we describe attempts to perform such an evaluation using a model-simulated 3-D CO 2 field as an independent reference.We constructed this reference field by running forward with NIES-TM a posteriori fluxes estimated from the GV data only (mentioned earlier in this section).The quality of this GV-based global CO 2 field was examined with the TCCON references.Figure 10 shows the monthly time series data collected at five TCCON sites (Ny Alesund, Norway; Białystok, Poland; Park Falls, USA; Tsukuba, Japan; and Wollongong, Australia; described in: Wunch et al., 2011a;Washenfelder et al., 2006;Messerschmidt et al., 2010Messerschmidt et al., , 2011;;Deutscher et al., 2010;Ohyama et al., 2009, respectively) and the corresponding forward simulation results.As summarized in Table 1, the difference between TCCON (version GGG 2009) and GV-based model predictions are mostly within the range of observational uncertainties.The addition of GOSAT retrievals decreases the differences at TCCON sites.The effects of differences in a priori concentration vertical profiles and column averaging kernels between TCCON and GOSAT retrievals were not taken into account, as their influence (estimated to be on the order of 0.1 ppm by Reuter et al., 2011) appear to be minor.Figure 11 compares the distribution of the 5 • × 5 • GOSAT X CO 2 values with the reference field (monthly-mean GOSAT X CO 2 minus the corresponding reference X CO 2 concentrations).Positive or negative deviations from the GV-based CO 2 reference greater than 2 ppm, which is larger than the range of TCCON uncertainty and modelobservation misfits, are found in lower South America (August and May), equatorial Africa (November and February), and central Asia (August and May).
When compared with the distributional patterns of X CO 2 values shown in Fig. 7, the locations of enhanced devia-tions from the reference coincide well with regions of higher UR values, which are indicative of greater involvement of GOSAT X CO 2 retrievals in the flux estimation.From the limited view point of this particular GV-based CO 2 concentration reference, we note that the changes in fluxes observed in Fig. 9 are likely induced by the GOSAT X CO 2 retrievals.This point is, of course, strictly subject to changes made in the current minimum observational uncertainty settings, which can increase the competition between GV data and GOSAT retrievals in constraining fluxes, or to improvements made in constructing the modeled CO 2 reference field (e.g.augmenting the GV data used in estimating the surface-based fluxes by adding other available surface observations).
To illustrate the effects of adding GOSAT retrievals to the GV data on the seasonality of the estimated fluxes, we plotted time series of fluxes for June 2009 to May 2010 for southern Africa (regions 21-24) and boreal Eurasia (region 25-28) in Fig. 12.These regions are underconstrained by the surface data and the fluxes change with the addition of GOSAT retrievals.The fluxes for most regions and seasons remained within natural variability bounds, which also served as a priori constraints.Winter fluxes of northern boreal Eurasia appear variable and underconstrained as the seasonality was not imposed in the a priori flux uncertainties.The practical value of GOSAT retrievals to the inverse modeling of surface fluxes can be confirmed if the fluxes estimated from GV and GOSAT retrievals are more accurate than those from GV only.Direct validation with independent regional flux estimates, such as those being developed by Canadell et al. (2011), is an attractive option, as they combine a number of top-down and bottom-up regional CO 2 flux estimates based on flux tower observations, terrestrial ecosystem models, forest carbon inventories, and multiple inverse model outputs (e. g.Gloor et al, 2012, Dolman et al, 2012).However, those flux data for our analysis period are not yet available.Another approach, common in inverse modeling, is to examine statistics of differences between the a priori and a posteriori fluxes (Tarantola, 2005).In our case we evaluated changes in the differences between the a priori and a posteriori fluxes brought by the addition of GOSAT retrievals.In 12 land regions (05, 06, 09, 10, 12, 14, 15, 26, 29, 32, 33, and 39), the addition of GOSAT retrievals brought the a posteriori fluxes closer to the a priori values, thereby reducing the differences between them.However, the opposite was seen in 7 regions (16, 17, 18, 22, 25, 27, and 42).It is not always clear if the deviations from the a priori fluxes are in the right or wrong direction.For instance, in regions 25 and 27 (southern boreal Asia), the addition of GOSAT retrievals increases summer uptake, while evidence indicates that use of our a priori fluxes in forward concentration simulations results in underestimation of the amplitude of CO 2 seasonal cycle observed in Siberia (Saeki et al., 2013).
Fig. 13 shows the regional flux changes relative to the a priori fluxes for those regions and months whose URs were larger than 20 % (a total of 80 monthly regional flux values).The change in the fluxes is expressed as where m GV is the flux estimated from GV only, and m GV+GOSAT is the flux estimated from GV and GOSAT retrievals.A positive value of m 2 means that the deviation from the a priori value is reduced (i.e., the a posteriori flux becomes closer to the a priori), and a negative value indicates that the addition of GOSAT retrievals moves the estimated flux away from the a priori value.Small changes in fluxes (∼ 1 PgC yr −1 region −1 ) show that the effect of adding GOSAT retrievals is statistically similar in both positive and negative changes.In contrast, large positive m 2 values indicate that the addition of GOSAT retrievals acts to suppress many large and likely erroneous deviations from the a priori fluxes in the underconstrained regions for which estimations are made with GV data alone.
As discussed by Gurney et al. (2002Gurney et al. ( , 2004)), the problem of estimating fluxes for unconstrained regions is not only limited to large flux uncertainties; departures of the a posteriori fluxes from a priori values vary largely with different transport models used even if observation data are the same.GOSAT-induced reductions in differences between the a posteriori and a priori fluxes seen in underconstrained regions may confirm two points.First, a reduction in the estimated flux uncertainty corresponds to a reduction in the deviation of the a posteriori fluxes from the a priori values.Second, GOSAT retrievals on average do not contradict the a priori flux estimates; otherwise, the addition of GOSAT retrievals to GV data would increase deviations of the a posteriori flux estimates from the a priori, while claiming a reduction in flux uncertainty at the same time.
Annual global total flux can be used for evaluating the validity of the estimated fluxes.For this 1 yr period (June 2009 to May 2010), the values of the global total flux, estimated from the GV data only and both the GV and GOSAT retrievals turned out to be 4.75 GtC yr −1 (land/ocean: −2.09/−2.00GtC yr −1 ) and 5.18 GtC yr −1 (land/ocean: −1.19/−2.48GtC yr −1 ), respectively.The fossil fuel emissions for the period was 8.64 GtC yr −1 , and the a priori uptakes by the terrestrial biosphere and the ocean are: −0.70/−1.71GtC yr −1 , respectively.The difference in global total fluxes based on the two sets of observational data may come from regionally and temporally varying differences between GOSAT retrievals and model-simulated X CO 2 values based on the GV-only a posteriori flux data, as seen in Fig. 11.The characterization of the temporal and spatial distribution of CO 2 growth rate for this period (2009-2010) is complicated by high temperatures observed in the Northern Hemisphere during the first half of 2010 (Galarneau et al.,  2012), which were followed by a heat wave and a large forest fire in western Russia (Barriopedro et al., 2011;Shvidenko et al., 2011).Further investigations are needed regarding the spatiotemporal variability of biases in the current X CO 2 retrievals.

Summary and conclusions
We developed a global carbon cycle modeling system designed specifically for the analysis of GOSAT X CO 2 retrievals and the estimation of regional CO 2 fluxes and their seasonal and interannual variabilities.The components for the forward modeling, NIES-TM, ODIAC fossil fuel emission inventory, VISIT terrestrial biosphere model, and the oceanic pCO 2 data assimilation system, were optimized to reproduce the seasonal, interannual, and spatial variability of atmospheric CO 2 and surface CO 2 fluxes.To simulate X CO 2 with better accuracy, stratospheric transport in NIES-TM was simulated on an isentropic grid, and the model was tuned to reproduce the age of stratospheric air.Using large point source data and night light observations, we improved the spatial distribution of fossil fuel CO 2 emissions.We confirmed that the ODIAC high-resolution emissions correlate well with detailed bottom up inventories at the horizontal resolutions of 0.5 • to 1 • .VISIT was optimized by simultaneously fitting selected physiological parameters to observed atmospheric CO 2 seasonal cycle, net primary production, and biomass distribution.Oceanic surface CO 2 fluxes were simulated with a 4-D variational data assimilation system which is forced by reanalyzed ocean currents and uses available pCO 2 observations and global oceanic pCO 2 climatology as constraints.A recent version of GFED data was used to account for fire emissions.
The capability of an inverse modeling scheme that optimizes monthly mean fluxes for 64 regions (Patra et al., 2005b) was expanded to handle GOSAT X CO 2 retrievals.The scheme takes into account the vertical sensitivity of GOSAT measurement stored in column averaging kernel data and the effects of a priori CO 2 vertical profiles on modelpredicted X CO 2 values.Globally constant offsets of GV data and GOSAT retrievals were estimated as optimized parameters.Screening of GOSAT retrievals was implemented using an observation-adjusted 3-D CO 2 concentration climatology.To confirm the consistency between the simulated variations of the ground based and total column data, model-predicted X CO 2 concentrations based on the GV-only a posteriori flux data were compared with TCCON observations, and they were in good agreement.
The inverse modeling scheme was then applied to the estimation of surface CO 2 fluxes from ground-based and aircraft data in the GV data and GOSAT FTS SWIR Level 2 X CO 2 data (ver.02.00).The a posteriori fluxes based on both GV and GOSAT data were found to be close to those based on GV data only in regions that are well constrained by GV data.On the other hand, the fluxes estimated for remote regions away from the GV sampling networks showed considerable changes when GOSAT X CO 2 retrievals were added.Sizeable reductions in flux uncertainties were observed in those cases.
The a posteriori fluxes for most regions and seasons fell within the range of natural flux variability estimated with the terrestrial and oceanic component models.The analysis of the a posteriori regional flux data and annual mean flux suggests that more improvements in the quantification and correction of biases in GOSAT X CO 2 retrievals are necessary.Also, there is a need for a coordinated effort to check and validate the regional monthly CO 2 flux estimates.

Supplementary material related to this article is available online at
. The inverse Published by Copernicus Publications on behalf of the European Geosciences Union.S. Maksyutov et al.: Regional CO 2 flux estimates for 2009-2010

Fig. 2 .
Fig. 2. Top: June 2009 to May 2010 averaged air-to-sea CO 2 prior fluxes (gC m −2 day −1 ) used in the inversion.Bottom: global integral of air-to-sea CO 2 fluxes (PgC yr −1 ) and corresponding global mean data uncertainties used in the inversion.

Fig. 3 .
Fig. 3. Global distribution of the annual mean CO 2 emissions due to burning fossil fuels.

Fig. 4 .
Fig. 4. Boundaries of the 64 source regions adopted in this study.The numbers on the figure are the region IDs.Regions shaded with dark blue are not considered in the flux estimation.

Fig. 5 .
Fig. 5.The number of GOSAT Level 2 X CO 2 data records per each of 5 • × 5 • grid cells during the months of August 2009, November 2009, February 2010, and May 2010.Red circles indicate the locations of the GV measurement sites chosen for this study.

Fig. 6 .
Fig. 6.GOSAT X CO 2 retrievals in the form of input to our inverse modeling scheme (gridded to 5 • × 5 • cells and averaged on a monthly time scale).Cells with three or more retrievals per month are shown here.The bias was corrected by raising each X CO 2 retrieval by 1.20 ppm.Overlaid are GLOBALVIEW values (in circles; monthly means).Values for the months of August 2009 (summer in the Northern Hemisphere), November 2009 (fall), February 2010 (winter), and May 2010 (spring) are shown.
presents the UR values for August 2009, November 2009, February 2010, and May 2010.As indicated in Eq. (4), the value of UR is affected by three factors: (1) the uncertainty in the observations and a priori fluxes, given by C D and C M , respectively; (2) the sensitivity of observations to surface fluxes (determined by atmospheric transport and stored in G); and

Fig. 7 .
Fig. 7. Percent reduction in the uncertainty of monthly surface flux estimates, attained by adding the GOSAT X CO 2 retrievals to the GLOB-ALVIEW dataset.

Fig. 8 .
Fig. 8. Monthly fluxes (g C m −2 day −1 ) estimated for the 64 subcontinental regions using GV data and GOSAT X CO 2 retrievals.Results for the months of August 2009 (summer in the Northern Hemisphere), November 2009 (fall), February 2010 (winter), and May 2010 (spring are shown.The values presented here are the sum of a priori fluxes (terrestrial biosphere exchange or ocean exchange + anthropogenic emissions + forest fire emissions) and the correction to the a priori flux determined via the optimization.Note the different color-coded scales used for land and ocean regions.

Fig. 9 .
Fig. 9. Differences between the fluxes estimated from GV data only and those from GV and GOSAT X CO 2 retrievals.Note that different color-coded scales are used for land and ocean regions.

Fig. 12 .
Fig. 12.Time series of monthly fluxes (g C m −2 day −1 ) for June 2009 to May 2010, for quadrants (top to bottom: SW, SE, NW, NE) in the South Africa (left column), boreal Eurasia (right column) subcontinental regions.The graphs show a prior fluxes (green lines), fluxes estimated from GV data alone (red lines), and fluxes estimated from GV and GOSAT retrievals (blue lines).The error bars show flux uncertainties.The gray vertical bars represent percent reduction in the uncertainty (UR, Eq. 6; scale on right side of graphs).Figures for all 64 regions are available in the supplementary information.

Fig. 13 .
Fig. 13.Changes in flux deviation from the a prior value due to the addition of the GOSAT retrievalks.The value is expressed as m 2(introduced in Eq. 7) in (PgC region −1 yr −1 ) 2 .Only the fluxes associated with significant uncertainty reductions (UR > 20 %) were considered here.

Table 1 .
Root-mean-square differences (RMS difference) between TCCON and modeled concentrations (in ppm) over one year between June 2009 and May 2010.Also listed is the RMS of TCCON observation uncertainty (TCCON uncertainty in ppm).