Ability of the 4-D-Var analysis of the GOSAT BESD XCO2 retrievals to characterize atmospheric CO2 at large and synoptic scales

This study presents results from the European Centre for Medium-Range Weather Forecasts (ECMWF) carbon dioxide (CO2) analysis system where the atmospheric CO2 is controlled through the assimilation of column-averaged dry-air mole fractions of CO2 (XCO2) from the Greenhouse gases Observing Satellite (GOSAT). The analysis is compared to a free-run simulation (without assimilation of XCO2), and they are both evaluated against XCO2 data from the Total Carbon Column Observing Network (TCCON). We show that the assimilation of the GOSAT XCO2 product from the Bremen Optimal Estimation Differential Optical Absorption Spectroscopy (BESD) algorithm during the year 2013 provides XCO2 fields with an improved mean absolute error of 0.6 parts per million (ppm) and an improved station-to-station bias deviation of 0.7 ppm compared to the free run (1.1 and 1.4 ppm, respectively) and an improved estimated precision of 1 ppm compared to the GOSAT BESD data (3.3 ppm). We also show that the analysis has skill for synoptic situations in the vicinity of frontal systems, where the GOSAT retrievals are sparse due to cloud contamination. We finally computed the 10-day forecast from each analysis at 00:00 UTC, and we demonstrate that the CO2 forecast shows synoptic skill for the largest-scale weather patterns (of the order of 1000 km) even up to day 5 compared to its own analysis. Abstract. This study presents results from the European Centre for Medium-Range Weather Forecasts (ECMWF) carbon dioxide (CO 2 ) analysis system where the atmospheric CO 2 is controlled through the assimilation of column-averaged dry-air mole fractions of CO 2 (XCO 2 ) from the Greenhouse gases Observing Satellite (GOSAT). The analysis is compared to a free-run simulation (without assimilation of XCO 2 ), and they are both evaluated against XCO 2 data from the Total Carbon Column Observing Network (TC-CON). We show that the assimilation of the GOSAT XCO 2 product from the Bremen Optimal Estimation Differential Optical Absorption Spectroscopy (BESD) algorithm during the year 2013 provides XCO 2 ﬁelds with an improved mean absolute error of 0.6 parts per million (ppm) and an improved station-to-station bias deviation of 0.7 ppm compared to the free run (1.1 and 1.4 ppm, respectively) and an improved estimated precision of 1 ppm compared to the GOSAT BESD data (3.3 ppm). We also show that the analysis has skill for synoptic situations in the vicinity of frontal systems, where the GOSAT retrievals are sparse due to cloud contamination. We ﬁnally computed the 10-day forecast from each analysis at 00:00 UTC, and we demonstrate that the CO 2 forecast shows synoptic skill for the largest-scale weather patterns (of the order of 1000 km) even up to day 5 compared to its own analysis.

Rather than using the relatively sparse network of the surface air-sample measurements, here we explore the measurements from satellite sounders in order to have a more global picture of the atmospheric CO 2 . To extract information on the CO 2 content in the atmosphere, passive atmospheric remote sounders measure in the thermal infrared (TIR) or in the near infrared/short-wave infrared (NIR/SWIR).
The Atmospheric Infrared Sounder (AIRS), measuring in the TIR, detects thermal radiation emitted by the Earth's surface and the atmosphere (Chédin et al., 2003). The assimilation of the AIRS observed radiances was developed by Engelen et al. (2009) at the European Centre for Medium-Range Weather Forecasts (ECMWF) using a four-dimensional variational (4-D-Var) data assimilation scheme. Their results showed the potential of data assimilation to constrain atmospheric CO 2 . They also showed the limitations of the assimilation of AIRS radiances, in particular due to the vertical sensitivity of the sounder. Due to the low thermal contrast between the Earth's surface and the air masses above, AIRS measurements have limited or no sensitivity to the lower troposphere and higher sensitivity to the middle atmosphere. Because the signals of the CO 2 surface sources and sinks are the largest in the near-surface and lower troposphere than in the middle atmosphere, AIRS measurements were not able to capture these signals.
In contrast, column-averaged dry-air mole fractions of CO 2 (or XCO 2 ) with a high near-surface sensitivity are retrieved from NIR/SWIR measurements based on scattered and back-scattered solar radiation; however, the NIR/SWIR measurements also have their limitations. They need sunlight and are therefore limited to daytime observations. Sufficiently cloud-free conditions and a low aerosol optical depth are also needed for accurate XCO 2 retrievals.
The aim of this study is to document the assimilation of XCO 2 products from NIR/SWIR measurements in order to constrain atmospheric CO 2 and to document how the assimilation impacts the simulated atmospheric CO 2 concentration. For that purpose, we assimilated the XCO 2 products derived from the NIR/SWIR spectra of the Greenhouse gases Observing Satellite (GOSAT; Kuze et. al, 2009). The assimilation system is based on the ECMWF system of Engelen et al. (2009), which has lately evolved for CAMS in order to assimilate retrieved products instead of observed radiances ).
The assimilation system provides an analysis of the atmospheric CO 2 concentration that is then integrated in time using a forecast model. The CO 2 forecast model used in this study is documented by Agustí-Panareda et al. (2014). In this model, the production and loss of CO 2 at the surface is based on surface fluxes that are partially prescribed and partially modelled. These CO 2 surface fluxes are not directly constrained by observations and they may deviate from reality. The accumulation of surface fluxes errors then leads to biases in the atmospheric CO 2 . On the other hand, the strength of the CO 2 forecast model is its ability to provide a realistic CO 2 synoptic variability. The first objective of this study is to determine the quality of the XCO 2 fields resulting from the assimilation of GOSAT XCO 2 data with a CO 2 forecast model where the CO 2 surface fluxes are not constrained.
The atmospheric CO 2 synoptic variability on a regional scale is related to the passage of frontal systems (Wang et al., 2007). These events are difficult to capture with the GOSAT measurements as the availability of the data is limited due to cloud contamination. Therefore, the second objective of this study is to document whether the assimilation helps improve the simulation of atmospheric CO 2 for synoptic events despite the lack of measurements nearby frontal systems.
Within CAMS, ECMWF is providing a CO 2 analysis based on the assimilation of GOSAT XCO 2 data with a delay of 5 days behind real time. A 10-day forecast is then issued from the analysis in order to provide the atmospheric CO 2 field in real time and for the next few days. The last objective of this study is to assess the quality of this forecast. The forecast quality as a function of the lead time and the season is evaluated against the analysis. This paper is structured as follows. Section 2 introduces the data sets used in this study. Section 3 describes our atmospheric CO 2 simulations with and without assimilation of the GOSAT XCO 2 data, and how we compared them with independent measurements. Sections 4 to 6 present the global evaluation of our simulations, a case study and the evaluation of the CO 2 forecast based on the analysis. Finally, Sect. 7 presents our conclusions.

Data sets
In this study, we used two sets of data. The first one is the measurements from the GOSAT's Fourier transform spectrometer and the XCO 2 product retrieved from these measurements by the University of Bremen (UoB) and is described in Sect. 2.1. The second one is the collection of measurements provided by the Total Carbon Column Observing Network (TCCON) and is described in Sect. 2.2.

GOSAT XCO 2
GOSAT is a joint effort between the Japanese Aerospace Exploration Agency (JAXA), the National Institute for Environmental Studies (NIES) and the Japanese Ministry of the Environment (MOE) as part of the Global Change Observation Mission (GCOM) programme of Japan. GOSAT was launched on 23 January 2009 and carries the thermal and near-infrared sensor for carbon observations, which consists of a Fourier transform spectrometer (TANSO-FTS) and a cloud and aerosol imager (TANSO-CAI).
In this study, we used XCO 2 retrieved from TANSO-FTS measurements of the upwelling radiance at the top of the atmosphere by the Bremen Optimal Estimation DOAS (differential optical absorption spectroscopy) (BESD) algorithm of UoB. The BESD algorithm was initially developed to retrieve XCO 2 from nadir measurements of the SCanning Imaging Absorption spectroMeter for Atmospheric CHar-tographY (SCIAMACHY) remote sensing spectrometer on the ENVIronment SATellite (ENVISAT; Reuter et al., 2010Reuter et al., , 2011. The BESD algorithm has been modified to also retrieve XCO 2 from GOSAT measurements. A detailed description of the GOSAT BESD algorithm can be found in Heymann et al. (2015). In brief, the algorithm uses three fitting windows, the O 2 -A band (12 920-13 195 cm −1 ), a weak CO 2 absorption band (6170-6278 cm −1 ) and a strong CO 2 band (4804-4896 cm −1 ) from both the medium-and highgain (respectively M-gain and H-gain) GOSAT nadir modes. An optimal estimation-based inversion technique is used to derive the most likely atmospheric state from every individual GOSAT measurement using a priori knowledge. The BESD algorithm explicitly accounts for atmospheric scattering by clouds and aerosols, reducing potential systematic biases. The scattering information on cloud and aerosols is mainly obtained from the O 2 -A and strong CO 2 absorption bands.
We used an inhomogeneous GOSAT BESD XCO 2 data set in this study as the GOSAT BESD algorithm was still under development. This intermediate version of the GOSAT BESD XCO 2 data is referred to as MACC GOSAT BESD XCO 2 (MACC standing for Monitoring Atmospheric Composition and Climate, the precursor of CAMS). Nevertheless, from the beginning of 2014 onwards, we have been assimilating the current version of the GOSAT BESD data (v01.00.02; Heymann et al., 2015) in near-real time.
The TANSO-FTS detector has a circular field of view of 10.5 km when projected on the Earth's surface (at exact nadir). In 2013, it measured in a mode with three measurements across track, and the footprints were separated by ∼ 263 km across track and ∼ 283 km along track. GOSAT can also operate in target mode, resulting in a finer sampling distance. For these specific situations, we further thinned the observations on a 1 • × 1 • grid by removing all the observations but one chosen at random. This procedure avoids having several measurements in the same model grid cell during the assimilation. This thinning, plus the characteristics of the instrument (measurement only during sunlit periods) and the processing of the level-2 data procedure (retrievals for clear-sky conditions and only over land), reduces the number of GOSAT XCO 2 data to about 100 per day. The assimilation window of 12 h means that about 50 GOSAT XCO 2 data points are assimilated during each time window.
The geographic distribution of these data is dependent on the season and the atmospheric conditions as illustrated by Fig. 1. For example, in July 2013, GOSAT BESD data are available up to 75 • N, and in October 2013 they are available only up to 60 • N. The reason for this is the solar geometry and the filtering of measurements under high solar zenith angle (SZA) conditions, where XCO 2 is more challenging to retrieve as the impact of atmospheric scattering becomes larger compared to low-SZA conditions. Other data gaps are due to the strict cloud filtering and other types of filtering, like those based on the quality of the spectral fits, on scattering parameters, on the meteorological state, and on the measurement geometry.
The MACC GOSAT BESD XCO 2 data sets have been bias-corrected using the TCCON data. As this data set is delivered in near-real time and the TCCON data are delivered with a delay of few months, it was not possible to directly compare the two data sets. Instead, the TCCON data from the previous year were used and they were corrected assuming a 2 parts per million (ppm) global atmospheric growth of CO 2 . A global offset was then computed and applied to the MACC GOSAT BESD XCO 2 based on the comparison between this data set and the corrected TCCON data set of the previous year. Moreover, with this procedure the TCCON data used in this study (same year as for the MACC GOSAT BESD XCO 2 data set) can be considered as independent data.
For the assimilation, the observation error covariances have to be specified. In this study, we assumed that the observation errors are not correlated in space and time. For the standard deviation of the observation error, we used the uncertainty of the BESD XCO 2 product provided together with the data. The BESD XCO 2 uncertainty product accounts for the various sources of uncertainty of the retrieval process. It varies in time and space around an average value of 2 ppm. We furthermore established that the specified observation error based on the XCO 2 uncertainty globally matches the expected observation error using diagnostics posterior to the analysis (not shown).

TCCON XCO 2
The TCCON is a network of ground-based Fourier transform spectrometers recording direct solar spectra in the near infrared spectral region (http://tccon.ornl.gov/). The columnaveraged dry-air mole fractions of CO 2 are retrieved from these spectra together with other chemical components of the atmosphere (Wunch et al., 2011a). In 2014, the version GGG2014 of the TCCON data was released. The errors on the retrieved XCO 2 are documented to be below 0.25 % (∼ 1 ppm) until the solar zenith angles are larger than 82 • .
When we downloaded the GGG2014 data in November 2015, 20 TCCON stations were providing data within the time period we are interested in (year 2013). Not all the stations were used in this study. First we removed JPL 2011 (USA), Pasadena/Caltech (USA) and Tsukuba (Japan), as they are not background stations and are associated with significant representativity errors. We also removed Edwards (USA). This station started to retrieve data from the middle of the year 2013, and we assumed that this was not long enough to provide information on the seasonal variation of the error in our simulations. Additionally, we removed Eureka (Canada) from the list of stations as the site was pro-  (Table 1). Orléans (France) had a specific treatment compared to the other stations. The averaging kernels were not specified in the GGG2014 release. Therefore, we decided to use the same information as for Lamont (USA) as advised in the previous release of the TCCON data (version GGG2012).

Experimental setup
We ran two model simulations for the year 2013. The first is similar to the operational CAMS CO 2 forecast  and is referred to as the "free run". This simulation is used as the reference to assess the impact of the assimilation of the GOSAT BESD XCO 2 data. The second simulation is the analysis in which the GOSAT XCO 2 data are assimilated and is referred to as the "analysis". The configuration of both simulations is described in Sect. 3.1. The simulations were evaluated against each other and also against the TCCON data. Section 3.2 introduces the method-ology used in comparison of simulations and the TCCON data.

Model simulations
The global simulations of atmospheric CO 2 are performed within the numerical weather prediction (NWP) framework of the Integrated Forecasting System (IFS). The CO 2 mass mixing ratio is directly transported within IFS as a tracer and is affected by surface fluxes. The transport is computed online and is updated each 12 h, benefiting from the assimilation of all the operational observations within the IFS 4-D-Var assimilation system. The terrestrial biogenic carbon fluxes are also computed online by the carbon module of the land surface model (Carbon-TESSEL or CTESSEL; Boussetta et al., 2013), while other prescribed fluxes are read from CO 2 surface fluxes inventories (see Agustí-Panareda et al., 2014, for more details).
The ability to assimilate retrieval products from GOSAT was included in IFS and is detailed in Massart et al. (2014) for the assimilation of methane data. The system used in this study is similar to the one of Massart et al. (2014) and is based on fixed background errors derived from the National Meteorological Center (NMC) method (Parrish and Derber, 1992). The standard deviation of the background error is constant for each model level and slowly increases from the upper troposphere to the lower troposphere with values from about 1 to about 5 ppm, and then rapidly increases to reach a value of about 40 ppm at the surface. The correlation of the background errors varies over the whole domain and vertically with a representative length scale of about 250 km. The system does not account for the spatial or temporal correlation between the errors of the observations. We chose in this study to have a horizontal resolution of TL255 on a reduced Gaussian grid (∼ 80 km×80 km), and 60 vertical levels from the surface up to 0.1 hPa. This resolution is sufficient for resolving the large-and synoptic-scale horizontal structures (∼ 1000 km) of the atmospheric CO 2 fields.

Comparison with TCCON
To evaluate the quality of the model simulations (free run and analysis), we have extensively used the TCCON data in this study. The comparison is performed in the TCCON space using the TCCON a priori and averaging kernel information (see Appendix A for more details). In order to have a decomposition of the errors of the model column-averaged CO 2 against the TCCON measurement, we computed for each TCCON station k for k ∈ [1, N ], the mean difference (or bias) δ k and the standard deviation of the difference (or scatter) σ k over the M k times t i for i ∈ [1, M k ] when we have a TCCON observation for the station k. Ifĉ o k (t i ) for i ∈ [1, M k ] is the observed TCCON XCO 2 time series for the station k, and ifĉ k (t i ) for i ∈ [1, M k ] is the model equivalent time se-ries, then the bias δ k and scatter σ k are defined by Additionally, we computed the correlation coefficient r k be- Heymann et al. (2015), we also computed the model offset δ, the mean absolute error (MAE) , the station-to-station bias deviation σ and the model precision π for the N TCCON stations The statistics for the comparisons of the simulations against the TCCON data have some gaps in time due to gaps in the availability of the TCCON data. They are also valid only where the TCCON sites are located, i.e. 16 points distributed over the globe. To have a more global overview of the model bias and scatter against the TCCON data, we smoothed these statistics in time and space (see Appendix B for more details). In summary, for the bias we averaged all the model-measurement differences for each TCCON site using a 1-week time bin. We then fit the time evolution of the weekly bias with a function that combines a linear and a harmonic component for each station. The second step is an extrapolation in space. For each week, the weekly biases

Global evaluation of the analysis
In this section we first present the characteristics of the XCO 2 derived from the free-run simulation when compared to the TCCON data. Second, we present the impact of the assimilation of the MACC GOSAT BESD XCO 2 comparing the XCO 2 from the analysis against the XCO 2 from the free run. Then, we discuss whether the analysis represents an improvement compared to the free run in terms of statistics against the TCCON data. Finally, we discuss the merits of the analysis compared to the MACC GOSAT BESD data using the TCCON data as a reference.

Free-run simulation vs. TCCON
When compared with the TCCON data, the free-run simulation has a mean offset δ of −0.36 ppm and a mean absolute error of 1.08 ppm (Table 2). However, the individual station bias δ k spans a range from 2.3 ppm at Ascension Island (Saint Helena, Ascension and Tristan da Cunha) to −2.9 ppm at Białystok (Poland). The station-to-station bias deviation σ of the free-run simulation then has a value of 1.27 ppm.
The variations of the bias as well as the seasonal cycle of the bias are highlighted in the Hovmöller diagram displayed in Fig. 2a. First, it shows that the initial condition of the free run has a positive bias of about 2 ppm over the tropical region (region between 23 • S and 23 • N) when compared to the TCCON data. This bias is reduced during the spring and reappears the next summer. It reaches its highest values in autumn with more than 2 ppm. These results are slightly different from those of Agustí-Panareda et al. (2014), where the model bias was found to be more constant in the tropical region when comparing the background CO 2 in the marine boundary layer with the National Oceanic and Atmospheric Administration (NOAA) GLOBALVIEW-CO2. Here, the evaluation of the bias in the tropics is driven by the comparison with XCO2 measurements from the TC-CON station of Ascension Island. For this station, the values of the bias from July to September result from the interpolation process as no measurements were reported during this period (Fig. S1 of the Supplement).
In contrast to the situation at the tropics, the initial condition of the free run has a negative bias at northern mid- Table 2. Statistics of the XCO 2 difference between the simulations (free run and analysis) and the average hourly TCCON data (model-TCCON): bias (δ k , in ppm), scatter (σ k , in ppm) and correlation coefficient (r k ). Also shown are the mean, the mean absolute error (MAE) and the deviation of the stations bias (respectively δ, and σ , in ppm), the mean scatter (π , in ppm) and the mean r (last three rows). The second column (N ) is the number of data used for computing the statistics.  (Fig. 4a). The negative bias at these mid-latitudes is nevertheless confirmed by the comparison with other stations, like Karlsruhe (Germany) and Park Falls (USA), where we have some data at the beginning of the year (Fig. 4b and c). The negative bias at northern mid-latitudes remains high during the whole year, with an absolute value generally greater than 1 ppm at the end of spring, and in June and December. This can be explained by the fact that the model does not release enough CO 2 before and after the growing season, i.e. March to May and October to December, and by the fact that, in the model, the onset of the CO 2 sink associated with the growing season starts too early in the season (Agustí-Panareda et al., 2014). The precision π of the free run measured by the average scatter between the simulation and the TCCON data is 1.4 ppm (Table 2). Similarly to the bias, the scatter varies in time and space as highlighted by the Hovmöller diagram of the scatter (Fig. 3a). The scatter has its highest values of more than 1 ppm at the northern mid-latitudes during May-June-July. This increase in the scatter is driven by the behaviour of the free run at Sodankylä. There, the simulation has a larger variability than the measurements. For example, at the end of June, the simulation presents a decrease of about 7 ppm in 36 h, whereas the measurements show a decrease of about 4 ppm (Fig. 4a). Elsewhere, there is also an increase in the scatter between May and July which is during the Northern Hemisphere growing season. This increase could be explained by the difficulty for CTESSEL to model the terrestrial biogenic carbon fluxes during the growing season, which leads to higher variability in the simulated atmospheric CO 2 .

Analysis vs. free run
To assess the impact of the assimilation of the MACC GOSAT BESD XCO 2 , we compared the evolution of XCO 2 from the analysis with XCO 2 from the free simulation. Figure 5 presents the Hovmöller diagram (time vs. latitude) of this difference. It shows that the first region where the analysis impacts XCO 2 is the tropics. There, compared to the free run, the analysis continuously decreases XCO 2 by up to 1 ppm in June and by more than 2 ppm from September to December. The assimilation of the GOSAT data consequently causes an improvement as the free run has a positive bias in this region in autumn compared to the TCCON data.
The analysis also decreases XCO 2 over the southern extra tropics (region between 23 and 66 • S) when compared to the free run (Fig. 5). The decrease extends to the southern high latitudes (≥ 66 • S) even when no GOSAT data were assimilated in this region. This decrease results mainly from the transport of CO 2 from the equatorial region and southern mid-latitudes towards southern high latitudes. Unfortunately, there are no independent XCO 2 data available at southern high latitudes to assess the merits of the analysis there.
Despite the fact that some GOSAT data are assimilated in the northern mid-latitudes during the first months of the simulation, the analysis only starts to differ significantly from the free run from March onwards. In this region, north of 30 • N, the analysis has higher values of XCO 2 than the free run, with a difference of more than 2 ppm during the northern summer. Again, the assimilation of the GOSAT data improves the simulated XCO 2 as the free run shows a strong negative bias there. Similar to the behaviour discussed for the southern high latitudes, the change in the CO 2 concentration at northern mid-latitudes is transported northward to higher latitudes. There is, nevertheless, a difference between the two hemispheres. For the Northern Hemisphere we have more data at high latitudes, especially during the summer, when the northernmost GOSAT measurements' cover goes up to 80 • N.

Analysis vs. TCCON data
When compared with the TCCON data, the GOSAT BESD XCO 2 analysis has an offset δ of −0.34 ppm and a mean absolute error of 0.57 ppm ( Table 2). The offset is similar to that of the free run (−0.36 ppm), but the mean absolute er-ror is improved (1.08 ppm for the free run). The individual station bias is moreover more constant in time for the analysis compared to the free run. For example, the trend of the free-run bias is 2.08 ppm yr −1 for Lauder (New Zealand) (Table S1 of the Supplement), and it improves to 0.47 ppm yr −1 for the analysis (Table S2 and Fig. 4c).
By increasing XCO 2 in the northern mid-latitudes as discussed before, the analysis considerably reduces the bias. A residual seasonal cycle in the bias is still present, with values usually in the range of 0 to 3 ppm (Fig. 2b). This could be explained by the fact that we correct the atmospheric state of CO 2 and not the CO 2 fluxes. During the seasons when the CO 2 fluxes are the main driver of the atmospheric CO 2 , the optimization of the atmospheric state only may not be enough.
The analysis has a more constant bias in time than the free run. It is also more accurate in space, with a station-to-station bias deviation σ that is largely reduced compared to the free run with a value of 0.61 ppm against 1.27 ppm ( Table 2). The assimilation of the MACC GOSAT BESD XCO 2 thus helps to significantly improve the accuracy of the model. The assimilation also helps improve the precision π , with the mean scatter improved by 15 %, reduced to a value of 1.22 ppm. The scatter of the analysis is reduced for all TCCON stations compared to the free run except for Garmisch (Germany), Zealand, between 1 January and 31 December 2013. For each station, the top panel presents the daily averaged data from TCCON (black dots), the daily averaged data from GOSAT co-located in time and space with the station (yellow squares), the simulated XCO 2 (solid lines) and the daily averaged simulated XCO 2 in the observation space (coloured dots). The bottom panel presents the weekly averaged bias of the simulated XCO 2 against the TCCON data (coloured dots) and the smoothed bias (solid lines). Blue represents the free run, while red is for the analysis.
where the scatter remains essentially unchanged. The Hovmöller diagram of the scatter shows that the main reduction is in the northern high latitudes in May (Fig. 3). In particular, the analysis shows less spurious variability than the free run at Sodankylä (Fig. 4a).

Analysis vs. MACC GOSAT BESD data
The analysis is much more accurate and more precise than the free run when compared to the TCCON data. The analysis also fills the gaps in time and space of the MACC GOSAT BESD data. In this section, we evaluate the analysis against the MACC GOSAT BESD data once more using the TCCON data as a reference.
The MACC GOSAT BESD data were compared to the TC-CON data using a geolocation criterion of 5 • in space and a time window of ±2 h. Before computing the difference between each GOSAT-TCCON pair, following Dils et. al (2014), we added a correction to the GOSAT-retrieved value in order to account for the use of different a priori CO 2 profiles in the two products. Moreover, we only kept the stations where more than 30 GOSAT-TCCON pairs were found in order to have more robust statistical results. This proce-dure removes Izaña (Spain), Ascension Island, Réunion Island (France) and Lauder from the list of the used TCCON stations in the comparison and reduces the number of stations to 12 (Table 3).
For each GOSAT-TCCON pair, we extracted the CO 2 profile from the analysis at the same location and time as the GOSAT measurement before computing the difference between the model and the TCCON data. In this way, we have a fair comparison between the analysis and the MACC GOSAT BESD data with respect to the TCCON data.
The resulting subset of the analysis minus TCCON differences has a different offset than the full data set but a similar mean absolute error, station-to-station bias deviation and precision (Tables 2 and 3). The difference in the offset is mainly due to a difference in the sampling between the subset and the full data set over the Northern Hemisphere. Due to few or no pairs occurring in spring for the subset, the sampling misses the negative bias of the analysis there. Missing the negative bias of the analysis results in an increased offset. In that respect, the mean absolute error is less sensitive to the used data set (subset or full data set).
The analysis has a lower mean absolute error than the one from the MACC GOSAT BESD data (0.65 ppm vs. Table 3. Statistics of the XCO 2 differences between the MACC GOSAT BESD data set and the average hourly TCCON data (left block, GOSAT-TCCON) or the analysis and the average hourly TCCON data (right block, model-TCCON): bias (δ k , in ppm), scatter (σ k , in ppm) and correlation coefficient (r k ). The analysis has been sampled similarly to the GOSAT data set in time and space. Also shown are the mean, the mean absolute error (MAE) and the deviation of the stations bias, the mean scatter (all in ppm) and the mean r (last three rows). The second column (N ) is the number of data points used for computing the statistics.  1 ppm, Table 3), a station-to-station bias deviation σ almost half of the one from GOSAT data (0.7 ppm vs. 1.3 ppm) and has an improved precision π (1 ppm vs. 3.3 ppm). The mean correlation coefficient is also higher in the analysis than in the satellite data with a value of 0.8 compared to 0.5. The statistics of the MACC GOSAT BESD data found here are different than those of Heymann et al. (2015), who used a more recent version of the GOSAT BESD product. With the successive improvements in the BESD algorithm, the latest version has a station-to-station bias deviation of ∼ 0.4 ppm and a precision of ∼ 2 ppm. The better precision (lower value of π) and the lower value of the mean absolute error and station-to-station bias deviation σ of the analysis compared to the MACC GOSAT BESD data set shows that the analysis is capable of smoothing the scatter of the satellite data. Moreover, the analysis is able to fill the gaps of the satellite data in time and space.

Case study of a cold front over Park Falls
The CO 2 concentration could be strongly affected by frontal systems. As an illustration, such a situation occurred at the end of May 2013, close to the TCCON station of Park Falls, Wisconsin, USA, when a cold front came from the northwest. On 31 May, the XCO 2 dropped from 398.62 ppm at 08:15 LT (local time) to 395.97 ppm at 12:53 LT (Fig. 6, top panel). This sudden decrease of 2.65 ppm in less than 5 h occurs after the arrival of a cold front, which is associated with a decrease in the surface pressure and a decrease in the temperature at 500 hPa (Fig. 6, lower panel).
The free run is able to capture the sudden decrease in XCO 2 , highlighting the skill of the model for such a situation (Fig. 6, upper panel). The flow during this period is mainly a descent of cold air from Canada towards the midwestern and eastern US. This cold air mass is depleted in CO 2 relative to the background (Figs. 7e and f). When it moves towards Park Falls, it results in decreasing XCO 2 as observed and simulated, but the decrease in the free run is too strong by 2 to 3 ppm compared to the measurements.
We investigated whether the assimilation of the GOSAT data helps improve the simulated evolution of the CO 2 concentration for such situations even if the number of BESD GOSAT data is limited in the vicinity of a frontal system due to the strict cloud filtering. Frontal systems are associated with clouds formed when moist air between the cold and warm fronts is lifted.
On 30 May, we have a few GOSAT measurements over the north and northeast region of North America (Fig. 7a). These measurements have the effect of increasing the XCO 2 in this region ( Fig. 7b-d). The cold air mass is then richer in CO 2 in the analysis compared to the free run, and when it moves towards Park Falls, the decrease is weaker and closer to the observed decrease. The assimilation of the GOSAT data helps improve the simulation by correcting the large-scale structure upstream and by improving the large-scale atmospheric XCO 2 horizontal gradient.
The XCO 2 decrease continues the next day on 1 June in both simulations as the cold front continued its descent. Unfortunately, likely due to the presence of clouds, no TCCON measurements are available during this period to corroborate the simulated XCO 2 decrease.

Forecast based on the analysis
Within CAMS, we are receiving the GOSAT BESD data for a given day with a delay of 5 days behind real time. The analysis for this day is run as soon as the data are received. A 10day forecast is then subsequently run based on the resulting analysis.
In this section, we aim to evaluate the forecast as a function of its lead time by comparing the forecast to the analysis valid for the same time. This comparison informs us about how long the information provided by the analysis remains in the forecast. Assuming perfect transport and perfect surface fluxes, the analysis and the forecast (valid for the same time) should be similar given that the analysis accurately corrects the atmospheric concentration of CO 2 . In practice, the differences observed between the analysis and the forecast could come from either the transport, the surface fluxes or the analysis.
To compare a forecast with the analysis valid for the same time, we computed the anomaly correlation coefficient (ACC) for XCO 2 (see Appendix C for more details). The ACC can be regarded as a skill score relative to the climatology: the higher the ACC, the better the forecast. In the framework of NWP, an ACC reaching 50 % corresponds to forecasts for which the error is the same as for a forecast based on a climatological average. An ACC of about 80 % indicates valuable skill in forecasting large-scale synoptic patterns.
We computed the ACC for each month individually as we know that the surface fluxes, drivers of the difference between the forecast and the analysis, have a strong seasonal cycle. We also computed it for different domains (globe, tropics and mid-to high latitudes) and for several forecast lead times, from 12 h up to 10 days. We found that the ACC is globally more than 90 % for day 3 and almost always more than 85 % for day 5 for each single month (Fig. 8a). This means that the forecast for today based on the analysis of 5 days ago shows the same large-scale synoptic XCO 2 patterns as the analysis. The information of the analysis therefore lasts long enough in the forecast to provide a good quality 5-day forecast for today (compared to the analysis). The information lasts longer in the tropics than in the Northern Hemisphere and slightly longer in the Northern Hemisphere than in the Southern Hemisphere ( Fig. 8b to d). This difference between the two hemispheres may reflect the fact that the CO 2 variability is much weaker in the Southern Hemisphere.
For forecasts longer than 5 days, globally, there are two particular months for which the ACC decreases faster than the others, i.e. July and December. For example, for these two months the ACC at day 5 is similar to the ACC at day 10 for October. This means that for July and December, the medium-range XCO 2 forecast (between 5 and 10 days) should be used more carefully. For July, the drop in skill occurs mainly over the Northern Hemisphere. The main reason is that the CO 2 fluxes are an even more important driver of the CO 2 concentration than the initial CO 2 concentration for this month. To better understand the impact of the surface fluxes, let us assume that in July we have too little release or, similarly, too much uptake of CO 2 in the atmosphere in the model over the Northern Hemisphere (as confirmed by Fig. 2a). This induces a negative bias of the CO 2 surface fluxes in the model. In the meantime, the analysis increases the CO 2 concentration helped by the GOSAT BESD data (Fig. 5). However, the next 12 h short-term forecast (used as the background for the next analysis) will not increase the CO 2 concentration enough due to the negative bias of the CO 2 fluxes. This opposition between the analysis and the short-term forecast explains the reduction in skill during the periods when the surfaces fluxes are the most important driver of the CO 2 concentration in the atmosphere. The global drop in skill for December is not directly related to a particular region as for July. It is nonetheless the second worst month for the tropics (after January) and the third worst for the Northern Hemisphere (together with September). Over the tropics during the winter, the reduction in skill is due to the opposite effect as for July over the Northern Hemisphere: the CO 2 fluxes are important and there is a positive bias in the fluxes (too much release or too little uptake of CO 2 in the atmosphere) in the model. For these situations when the CO 2 fluxes are the main driver of the atmospheric CO 2 , the only solution to improve the skill would be to optimize the CO 2 fluxes together with the CO 2 initial conditions.

Conclusions
The Copernicus Atmosphere Monitoring Service (CAMS) greenhouse gases data assimilation within the numerical weather prediction (NWP) framework of the Integrated Forecasting System (IFS) is designed to correct the atmospheric concentration of CO 2 instead of the surface fluxes in order to constrain the atmospheric CO 2 . This requires the use of a short assimilation window so as to neglect the model errors of the short-term forecast (lasting the length of the assimilation window). In the case of atmospheric CO 2 , model errors are related to potentially inaccurate surface fluxes or transport.
This article demonstrates the benefit of the assimilation of XCO 2 data derived from the Greenhouse gases Observing Satellite (GOSAT) by intermediate versions of the Bremen Optimal Estimation DOAS (BESD) algorithm of the University of Bremen (UoB). The assimilation of the GOSAT BESD XCO 2 provides a CO 2 analysis that was compared to a free-run forecast where the CO 2 concentration is not constrained by any CO 2 observation. The comparison was 1 year long (year 2013) and both simulations (analysis and free run) were evaluated against measurements from the Total Carbon Column Observing Network (TCCON). We showed that the free run has a negative bias at northern mid-latitudes and a large positive bias in the tropical region with strong seasonal variations in both regions. These results are consistent with the biases documented by Agustí-Panareda et al. (2014) and mainly associated with biogenic fluxes.
The analysis significantly reduces these biases without completely removing them, with a remaining mean offset of −0.34 ppm and a mean absolute error of 0.57 ppm compared to the TCCON data. However, the accuracy estimated with the station-to-station bias deviation is 0.61 ppm. This represents a large improvement compared to the free run, for which the accuracy is 1.27 ppm. The precision of the analysis estimated with the mean scatter is 1.22 ppm, slightly better than for the free run with a value of 1.43 ppm.
The analysis produced in this paper was compared to the assimilated MACC GOSAT BESD data using TCCON data as a reference. This comparison showed that the analysis has a lower station-to-station bias deviation than the assimilated data (0.7 ppm compared to 1.3 ppm). The precision is much better for the analysis, with a scatter of 1 ppm, while the assimilated data have a scatter of 3.3 ppm. The precision of the analysis is also better than the documented precision of other GOSAT XCO 2 products. The precision of the NIES product extracted from Yoshida et. al (2013) is 1.8 ppm. The precision of the University of Leicester product and of the SRON Netherlands Institute for Space Research product is respectively 2.5 and 2.37 ppm (Dils et. al, 2014). The CO 2 analysis is consequently an alternative to the standard XCO 2 GOSAT products as it provides a lower or similar station-to-station bias deviation and a better-precision XCO 2 product compared to TCCON. Moreover, it has a uniform spatio-temporal resolution.
The pre-operational CAMS CO 2 analysis is similar to the analysis presented in this paper, having nevertheless a higher horizontal resolution (TL511 on a reduced Gaussian grid, ∼ 40 km×40 km), and a higher vertical resolution with 137 vertical levels. It currently assimilates the most recent version of the GOSAT BESD data presented by Heymann et al. (2015) in near-real time. These data have an improved bias deviation (∼ 0.4 ppm) and an improved precision (∼ 2 ppm) compared to those used in this study. The near-real-time CAMS CO 2 analysis should therefore have an improved station-to-station bias deviation and precision than the analysis presented in this paper.
We corrected the atmospheric concentration by only constraining the atmospheric concentration and not the surface fluxes. When and where the surface flux is a significant driver of the atmospheric concentration and if the assimilated data are not good enough or not numerous enough (in time and space), then constraining only atmospheric CO 2 does not compensate for the error in the surface flux. The next step is to further improve the carbon module CTESSEL in order to reduce the bias of the model. Another long-term solution would be to constrain the surface flux at the same time as the concentration.
One strength of the CO 2 model used in this study is its ability to represent CO 2 variations associated with synoptic weather systems (Agustí-Panareda et al., 2014). By correcting the large-scale XCO 2 patterns and removing part of the model bias, we showed with a case study that the analysis is able to better represent the CO 2 variations associated with these situations. The variations in the atmospheric reservoir of CO 2 are the result of changes in the surface fluxes to and from the atmosphere. If the characteristics of the analysis are found to be satisfactory in terms of bias and precision, the analysis could be included into a flux inversion system to infer surface fluxes.
The horizontal resolution of this study is half the horizontal resolution of the pre-operational analysis and the vertical resolution of the pre-operational analysis is also higher. One should expect an even better representation of the CO 2 variability in the pre-operational analysis. In the future, the horizontal resolution could be increased even further toward the ECMWF operational resolution of about 16 km × 16 km.
The quality of the analysis is considered to be sufficient to assess the quality of the forecast as a function its lead time. We showed that the forecast for day 3 and day 5, which will be the valid range for today's forecast, has an anomaly correlation coefficient of 90 and 85 %, respectively. This means that we are providing a CO 2 forecast with accurate synoptic features for today. With a good representation of the variability and a bias mostly under 1 ppm, the CAMS atmospheric CO 2 promises to become a useful product, for example, for planning a measurement campaign. It could also be used as the a priori in the satellite or TCCON retrieval algorithms or be used to evaluate the retrieval products from the Orbiting Carbon Observatory-2 (OCO-2, oco.jpl.nasa.gov).

Appendix A: Comparing the model against TCCON
For the comparison with the TCCON data, one has to account for the a priori information used in the retrieval that linksĉ o , the TCCON-retrieved XCO 2 to x t , the true (unknown) CO 2 profile (Wunch et al., 2011b), where x b is an a priori profile of CO 2 , a is a vector resulting from the product of the averaging kernel matrix with a drypressure weighting function vector (for the vertical integration), c b is the column-averaged mixing ratio computed from x b , and ε is the error in the retrieved column-averaged mixing ratio. This error includes the random and systematic errors in the measured signal and in the retrieval algorithm. To compare the model with the TCCON-retrieved value, we used the same a priori information, so that the model profile x is converted to a column-averaged mixing ratioĉ bŷ The comparison between the simulation and TCCON occurs in the observation space with the difference between the model column-averaged mixing ratioĉ of Eq. (A2) and the TCCON column-averaged mixing ratioĉ o of Eq. (A1), Let us define η = a T (x − x t ) as the model error in terms of the column-averaged mixing ratio. It accounts for numerous errors, for example the errors directly linked to the model processes like the transport, the errors in the surface fluxes, the representativity error and the error due to the assimilation of the GOSAT XCO 2 data for the analysis. The difference between the smooth model column-averaged mixing ratioĉ and the TCCON column-averaged mixing ratioĉ o is, therefore, the sum of the model error η and the error in the retrieved column-averaged mixing ratio ε.
To compute the model column-averaged mixing ratioĉ of Eq. (A2) equivalent to each TCCON measurement, we extracted the two model profiles that are closest to the measurement time and at the nearest grid point to the measurement. The two profiles are then interpolated in time in order to obtain the model profile at the same time as the measurement. Finally, we computed the column-averaged mixing ratio according to Eq. (A2).

Appendix B: Smoothing the statistics against TCCON
In order to have a more global view of the bias and the scatter of a simulation against the data from the TCCON network, we have developed and used a two-step algorithm. The first step consists in computing the statistics (bias and the standard deviation) for each week of 2013 and for each TCCON station when the data are available. The weekly statistics are then interpolated in time using a function described in the following section (Sect. B1). This allows one to fill in the gaps in time when no data are available. We therefore have a value for the bias at each station and for each week. For the second step, we compute a quadratic function of latitude that best fits the interpolated biases for each week (Sect. B2).

B1 Time smoothing
For each TCCON station k and for each week w l for l ∈ [1, 52], we compute the mean difference δ l k and the standard deviation of the difference σ l k between every TCCON observation during this week and the model equivalent value. The statistics are computed only when more than 10 TCCON measurements are available during the week. The averaged difference (or bias) is then interpolated in time t with the function b k (t) that combines a linear growth and a harmonic component, b k (t) = a k t + b k + α k sin t τ 1 + ϕ k + β k sin t τ 2 + ϕ k .
(B1) a k , b k , α k , β k and ϕ k are the parameters of the function b k (t) obtained by an optimization procedure that minimizes the distance between b k (t) and the series of δ l k for l ∈ [1, 52]. τ 1 is chosen to be 6 months and τ 2 3 months. The form of the function of Eq. (B1) thus gives a linear growing bias and allows seasonal variations. A similar function is used for the standard deviation.

B2 Spatial smoothing
The time smoothing allows us to fill in the gaps in the time series of the bias for each station, when for a given week we do not have any measurement to compare with. Following Bergamaschi et al. (2009), we then compute for each week w l the best fit of the interpolated biases with a quadratic function of latitudeb l , b l (φ) = a l φ 2 + b l φ + c l , where φ is the sine of the latitude. a l , b l and c l are obtained by an optimization procedure that minimizes the distance betweenb l and the weekly interpolated biases δ l k for k ∈ [1, N ]. A similar function is used for the standard deviation.

B3 Discussion
For some stations, the availability of the weekly differences is not uniform in time and the time smoothing of Eq. (B1) provides spurious values. We solved this issue by fixing the coefficient α k to a zero value (See Table S1).
With a root mean square error (RMSE) mostly under 0.7 ppm and a correlation mostly over 0.8, the smoothed bias matches well with the weekly bias (Table S1). The Hovmöller diagram (Fig. 2) can, thus, be considered as an accurate representation of the overall bias.
Compared to the bias, the fit between the time series of the weekly scatter and the regression is not as good for the scatter. The correlation coefficient is mostly between 0.5 and 0.7 (Table S1).
The Supplement related to this article is available online at doi:10.5194/acp-16-1653-2016-supplement.