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

Abstract

This study presents results from the European Centre for Medium-Range Weather Forecasts (ECMWF) carbon dioxide (${\mathrm{CO}}_{\mathrm{2}}$) analysis system where the atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ is controlled through the assimilation of column-averaged dry-air mole fractions of ${\mathrm{CO}}_{\mathrm{2}}$ (${\mathrm{XCO}}_{\mathrm{2}}$) from the Greenhouse gases Observing Satellite (GOSAT). The analysis is compared to a free-run simulation (without assimilation of ${\mathrm{XCO}}_{\mathrm{2}}$), and they are both evaluated against ${\mathrm{XCO}}_{\mathrm{2}}$ data from the Total Carbon Column Observing Network (TCCON). We show that the assimilation of the GOSAT${\mathrm{XCO}}_{\mathrm{2}}$ product from the Bremen Optimal Estimation Differential Optical Absorption Spectroscopy (BESD) algorithm during the year 2013 provides${\mathrm{XCO}}_{\mathrm{2}}$ fields with an improved mean absolute error of 0.6 parts per million ($\mathrm{ppm}$) and an improved station-to-station bias deviation of 0.7 $\mathrm{ppm}$ compared to the free run (1.1 and 1.4 $\mathrm{ppm}$, respectively) and an improved estimated precision of 1 $\mathrm{ppm}$ compared to the GOSAT BESD data (3.3 $\mathrm{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 $\mathrm{UTC}$, and we demonstrate that the ${\mathrm{CO}}_{\mathrm{2}}$ forecast shows synoptic skill for the largest-scale weather patterns (of the order of 1000 $\mathrm{km}$) even up to day 5 compared to its own analysis.

Introduction

Carbon in the atmosphere is present mostly in the form of carbon dioxide (${\mathrm{CO}}_{\mathrm{2}}$). Its amount is relatively small compared to the amount of carbon present in other reservoirs like the ocean . Being well mixed, atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ is nevertheless easier to monitor by measurements than other carbon reservoirs. To improve the monitoring of atmospheric ${\mathrm{CO}}_{\mathrm{2}}$, one can combine atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ measurements with a numerical model. This paper describes such a system, which has been developed for the Copernicus Atmosphere Monitoring Service (CAMS).

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 ${\mathrm{CO}}_{\mathrm{2}}$. To extract information on the ${\mathrm{CO}}_{\mathrm{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 . The assimilation of the AIRS observed radiances was developed by 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 ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$ (or ${\mathrm{XCO}}_{\mathrm{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 ${\mathrm{XCO}}_{\mathrm{2}}$ retrievals.

The aim of this study is to document the assimilation of ${\mathrm{XCO}}_{\mathrm{2}}$ products from NIR/SWIR measurements in order to constrain atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ and to document how the assimilation impacts the simulated atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ concentration. For that purpose, we assimilated the ${\mathrm{XCO}}_{\mathrm{2}}$ products derived from the NIR/SWIR spectra of the Greenhouse gases Observing Satellite GOSAT;. The assimilation system is based on the ECMWF system of , 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 ${\mathrm{CO}}_{\mathrm{2}}$ concentration that is then integrated in time using a forecast model. The ${\mathrm{CO}}_{\mathrm{2}}$ forecast model used in this study is documented by . In this model, the production and loss of ${\mathrm{CO}}_{\mathrm{2}}$ at the surface is based on surface fluxes that are partially prescribed and partially modelled. These ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$. On the other hand, the strength of the ${\mathrm{CO}}_{\mathrm{2}}$ forecast model is its ability to provide a realistic ${\mathrm{CO}}_{\mathrm{2}}$ synoptic variability. The first objective of this study is to determine the quality of the ${\mathrm{XCO}}_{\mathrm{2}}$ fields resulting from the assimilation of GOSAT ${\mathrm{XCO}}_{\mathrm{2}}$ data with a ${\mathrm{CO}}_{\mathrm{2}}$ forecast model where the ${\mathrm{CO}}_{\mathrm{2}}$ surface fluxes are not constrained.

The atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ synoptic variability on a regional scale is related to the passage of frontal systems . 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 ${\mathrm{CO}}_{\mathrm{2}}$ for synoptic events despite the lack of measurements nearby frontal systems.

Within CAMS, ECMWF is providing a ${\mathrm{CO}}_{\mathrm{2}}$ analysis based on the assimilation of GOSAT ${\mathrm{XCO}}_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{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  introduces the data sets used in this study. Section  describes our atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ simulations with and without assimilation of the GOSAT ${\mathrm{XCO}}_{\mathrm{2}}$ data, and how we compared them with independent measurements. Sections  to present the global evaluation of our simulations, a case study and the evaluation of the ${\mathrm{CO}}_{\mathrm{2}}$ forecast based on the analysis. Finally, Sect.  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 ${\mathrm{XCO}}_{\mathrm{2}}$ product retrieved from these measurements by the University of Bremen (UoB) and is described in Sect. . The second one is the collection of measurements provided by the Total Carbon Column Observing Network (TCCON) and is described in Sect. .

GOSAT ${\mathrm{XCO}}_{\mathrm{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 ${\mathrm{XCO}}_{\mathrm{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 ${\mathrm{XCO}}_{\mathrm{2}}$ from nadir measurements of the SCanning Imaging Absorption spectroMeter for Atmospheric CHartographY (SCIAMACHY) remote sensing spectrometer on the ENVIronment SATellite ENVISAT;. The BESD algorithm has been modified to also retrieve ${\mathrm{XCO}}_{\mathrm{2}}$ from GOSAT measurements. A detailed description of the GOSAT BESD algorithm can be found in . In brief, the algorithm uses three fitting windows, the ${\mathrm{O}}_{\mathrm{2}}$-A band (12 920–13 195 ${\mathrm{cm}}^{-\mathrm{1}}$), a weak ${\mathrm{CO}}_{\mathrm{2}}$ absorption band (6170–6278 ${\mathrm{cm}}^{-\mathrm{1}}$) and a strong ${\mathrm{CO}}_{\mathrm{2}}$ band (4804–4896 ${\mathrm{cm}}^{-\mathrm{1}}$) from both the medium- and high-gain (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 ${\mathrm{O}}_{\mathrm{2}}$-A and strong ${\mathrm{CO}}_{\mathrm{2}}$ absorption bands.

We used an inhomogeneous GOSAT BESD ${\mathrm{XCO}}_{\mathrm{2}}$ data set in this study as the GOSAT BESD algorithm was still under development. This intermediate version of the GOSAT BESD ${\mathrm{XCO}}_{\mathrm{2}}$ data is referred to as MACC GOSAT BESD ${\mathrm{XCO}}_{\mathrm{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; in near-real time.

The TANSO-FTS detector has a circular field of view of 10.5 $\mathrm{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 $\sim 263$ $\mathrm{km}$ across track and $\sim 283$ $\mathrm{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 ${\mathrm{1}}^{\circ }\times {\mathrm{1}}^{\circ }$ 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 ${\mathrm{XCO}}_{\mathrm{2}}$ data to about 100 per day. The assimilation window of 12 h means that about 50 GOSAT ${\mathrm{XCO}}_{\mathrm{2}}$ data points are assimilated during each time window.

Example of the distribution of the assimilated GOSAT BESD ${\mathrm{XCO}}_{\mathrm{2}}$ data: July 2013 (top panel, about 3400 retrievals) and October 2013 (bottom, about 1270 retrievals). The monthly data are here aggregated on a ${\mathrm{2}}^{\circ }\times {\mathrm{2}}^{\circ }$ grid and averaged. The blue/red represents the low/high averaged ${\mathrm{XCO}}_{\mathrm{2}}$ values in $\mathrm{ppm}$.

[Figure omitted. See PDF]

The geographic distribution of these data is dependent on the season and the atmospheric conditions as illustrated by Fig. . For example, in July 2013, GOSAT BESD data are available up to 75${}^{\circ }$ N, and in October 2013 they are available only up to 60${}^{\circ }$ N. The reason for this is the solar geometry and the filtering of measurements under high solar zenith angle (SZA) conditions, where ${\mathrm{XCO}}_{\mathrm{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 ${\mathrm{XCO}}_{\mathrm{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 ($\mathrm{ppm}$) global atmospheric growth of ${\mathrm{CO}}_{\mathrm{2}}$. A global offset was then computed and applied to the MACC GOSAT BESD ${\mathrm{XCO}}_{\mathrm{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 ${\mathrm{XCO}}_{\mathrm{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 ${\mathrm{XCO}}_{\mathrm{2}}$ product provided together with the data. The BESD ${\mathrm{XCO}}_{\mathrm{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 $\mathrm{ppm}$. We furthermore established that the specified observation error based on the ${\mathrm{XCO}}_{\mathrm{2}}$ uncertainty globally matches the expected observation error using diagnostics posterior to the analysis (not shown).

TCCON ${\mathrm{XCO}}_{\mathrm{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 column-averaged dry-air mole fractions of ${\mathrm{CO}}_{\mathrm{2}}$ are retrieved from these spectra together with other chemical components of the atmosphere . In 2014, the version GGG2014 of the TCCON data was released. The errors on the retrieved ${\mathrm{XCO}}_{\mathrm{2}}$ are documented to be below 0.25 % ($\sim \mathrm{1}$ $\mathrm{ppm}$) until the solar zenith angles are larger than 82${}^{\circ }$ (Wunch et al., 2015).

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 providing data during only 3 days in 2013. This selection of the TCCON stations left 16 stations for the study (Table ).

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).

List of the TCCON stations used, ordered by latitude from north to south.

Site	Lat	Long	Starting date	Reference
Sodankylä (sodankyla01)	67.37	26.63	6 Feb 2009
Białystok (bialystok01)	53.23	23.02	1 Mar 2009
Bremen (bremen01)	53.10	8.85	6 Jan 2005
Karlsruhe (karlsruhe01)	49.10	8.44	19 Apr 2010
Orléans (orleans01)	47.97	2.11	29 Aug 2009
Garmisch (garmisch01)	47.48	11.06	16 Jul 2007
Park Falls (parkfalls01)	45.94	$-$90.27	26 May 2004
Four Corners (fourcorners01)	36.80	$-$108.48	1 Mar 2011
Lamont (lamont01)	36.60	$-$97.49	6 Jul 2008
Saga (saga01)	33.24	130.29	28 Jul 2011
Izaña (izana01)	28.30	$-$16.48	18 May 2007
Ascension Island (ascension01)	$-$7.92	$-$14.33	22 May 2012
Darwin (darwin01)	$-$12.43	130.89	28 Aug 2005
Réunion Island (reunion01)	$-$20.90	55.49	6 Oct 2011
Wollongong (wollongong01)	$-$34.41	150.88	26 Jun 2008
Lauder 125HR (lauder02)	$-$45.05	169.68	2 Feb 2010

Experimental setup

We ran two model simulations for the year 2013. The first is similar to the operational CAMS ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{XCO}}_{\mathrm{2}}$ data. The second simulation is the analysis in which the GOSAT ${\mathrm{XCO}}_{\mathrm{2}}$ data are assimilated and is referred to as the “analysis”. The configuration of both simulations is described in Sect. . The simulations were evaluated against each other and also against the TCCON data. Section  introduces the methodology used in comparison of simulations and the TCCON data.

Model simulations

The global simulations of atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ are performed within the numerical weather prediction (NWP) framework of the Integrated Forecasting System (IFS). The ${\mathrm{CO}}_{\mathrm{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 $\mathrm{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;, while other prescribed fluxes are read from ${\mathrm{CO}}_{\mathrm{2}}$ surface fluxes inventories (see , for more details).

The ability to assimilate retrieval products from GOSAT was included in IFS and is detailed in for the assimilation of methane data. The system used in this study is similar to the one of and is based on fixed background errors derived from the National Meteorological Center (NMC) method . 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 $\mathrm{ppm}$, and then rapidly increases to reach a value of about 40 $\mathrm{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 $\mathrm{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 ($\sim 80\mathrm{km}\times 80\mathrm{km}$), and 60 vertical levels from the surface up to 0.1 $\mathrm{hPa}$. This resolution is sufficient for resolving the large- and synoptic-scale horizontal structures ($\sim 1000\mathrm{km}$) of the atmospheric ${\mathrm{CO}}_{\mathrm{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  for more details). In order to have a decomposition of the errors of the model column-averaged ${\mathrm{CO}}_{\mathrm{2}}$ against the TCCON measurement, we computed for each TCCON station $k$ for $k\in \left[\mathrm{1},,,N\right]$, the mean difference (or bias) ${\mathit{\delta }}_k$ and the standard deviation of the difference (or scatter) ${\mathit{\sigma }}_k$ over the $M_k$ times $t_i$ for $i\in \left[\mathrm{1},,,M_k\right]$ when we have a TCCON observation for the station $k$. If ${\widehat{c}}_k^o\left(t_i\right)$ for $i\in \left[\mathrm{1},,,M_k\right]$ is the observed TCCON ${\mathrm{XCO}}_{\mathrm{2}}$ time series for the station $k$, and if ${\widehat{c}}_k\left(t_i\right)$ for $i\in \left[\mathrm{1},,,M_k\right]$ is the model equivalent time series, then the bias ${\mathit{\delta }}_k$ and scatter ${\mathit{\sigma }}_k$ are defined by $$\begin{array}{cc}& {\mathit{\delta }}_k={\displaystyle \frac{\mathrm{1}}{M_k}}\sum _{i=\mathrm{1}}^{M_k}\left[{\widehat{c}}_k,\left(t_i\right),-,{\widehat{c}}_k^o,\left(t_i\right)\right],\\ & & {\mathit{\sigma }}_k=\sqrt{{\displaystyle \frac{\mathrm{1}}{M_k-\mathrm{1}}}\sum _{i=\mathrm{1}}^{M_k}{\left[{\widehat{c}}_k,\left(t_i\right),-,{\widehat{c}}_k^o,\left(t_i\right),-,{\mathit{\delta }}_k\right]}^{\mathrm{2}}}.\end{array}$$ Additionally, we computed the correlation coefficient $r_k$ between ${\widehat{c}}_k\left(t_i\right)$ and ${\widehat{c}}_k^o\left(t_i\right)$ for $i\in \left[\mathrm{1},,,M_k\right]$.

Hovmöller diagram (latitude vs. time) of the smoothed bias (in ppm, negative/positive in blue/red) of the simulated ${\mathrm{XCO}}_{\mathrm{2}}$ against the data of the TCCON network, from 1 January to 31 December 2013. (a) Free-run simulation. (b) Analysis.

[Figure omitted. See PDF]

Following , we also computed the model offset $\mathit{\delta }$, the mean absolute error (MAE) $\mathrm{\Delta }$, the station-to-station bias deviation $\mathit{\sigma }$ and the model precision $\mathit{\pi }$ for the $N$ TCCON stations $$\begin{array}{ccc}& \mathit{\delta }={\displaystyle \frac{\mathrm{1}}{N}}\sum _{k=\mathrm{1}}^N{\mathit{\delta }}_k,& \mathrm{\Delta }={\displaystyle \frac{\mathrm{1}}{N}}\sum _{k=\mathrm{1}}^N\left|{\mathit{\delta }}_k\right|,\\ & & \mathit{\sigma }=\sqrt{{\displaystyle \frac{\mathrm{1}}{N-\mathrm{1}}}\sum _{k=\mathrm{1}}^N{\left[{\mathit{\delta }}_k,-,\mathit{\delta }\right]}^{\mathrm{2}}},& \mathit{\pi }={\displaystyle \frac{\mathrm{1}}{N}}\sum _{k=\mathrm{1}}^N{\mathit{\sigma }}_k.\end{array}$$

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  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 of every station are extrapolated using a quadratic function of latitude. This results in a Hovmöller diagram of the bias as a function of time and latitude. A similar process is applied for the scatter (see Figs.  and ).

Statistics of the ${\mathrm{XCO}}_{\mathrm{2}}$ difference between the simulations (free run and analysis) and the average hourly TCCON data (model–TCCON): bias (${\mathit{\delta }}_k$, in $\mathrm{ppm}$), scatter (${\mathit{\sigma }}_k$, in $\mathrm{ppm}$) and correlation coefficient ($r_k$). Also shown are the mean, the mean absolute error (MAE) and the deviation of the stations bias (respectively $\mathit{\delta }$, $\mathrm{\Delta }$ and $\mathit{\sigma }$, in $\mathrm{ppm}$), the mean scatter ($\mathit{\pi }$, in $\mathrm{ppm}$) and the mean $r$ (last three rows). The second column ($N$) is the number of data used for computing the statistics.

		Free run				Analysis
Site	$N$	Bias	Scatter	$r$		Bias	Scatter	$r$
Sodankylä	20441	$-$1.59	1.35	0.91		$-$0.55	1.35	0.92
Białystok	16063	$-$2.68	1.96	0.81		$-$1.66	1.80	0.77
Bremen	4883	$-$1.62	1.52	0.79		$-$0.41	1.27	0.82
Karlsruhe	4201	$-$1.26	1.72	0.80		$-$0.25	1.54	0.82
Orléans	8444	$-$0.38	1.36	0.85		0.09	1.21	0.91
Garmisch	10371	$-$0.92	1.59	0.82		$-$0.29	1.62	0.80
Park Falls	27991	$-$1.69	2.06	0.81		$-$0.60	1.45	0.90
Four Corners	19872	0.69	1.76	0.58		0.57	1.43	0.74
Lamont	43731	$-$0.20	2.09	0.59		$-$0.04	1.35	0.80
Saga	10349	$-$1.19	1.61	0.75		$-$0.64	1.33	0.83
Izaña	4463	0.27	0.80	0.90		0.40	0.62	0.94
Ascension Island	7111	2.31	1.29	0.24		0.72	1.27	0.21
Darwin	29194	1.57	1.12	0.78		$-$0.02	1.04	0.79
Réunion Island	18880	0.56	0.73	0.76		$-$0.77	0.60	0.78
Wollongong	27562	0.30	1.05	0.71		$-$1.08	1.06	0.65
Lauder	53500	0.01	0.83	0.86		$-$0.97	0.59	0.85
Mean	16	$-$0.36	1.43	0.75		$-$0.34	1.22	0.78
MAE	16	1.08	–	–		0.57	–	–
Deviation	16	1.27	–	–		0.61	–	–

Same as Fig.  but for the standard deviation, with yellow/red for low/high values.

[Figure omitted. See PDF]

Time series of ${\mathrm{XCO}}_{\mathrm{2}}$ (in ppm) at (a) Sodankylä, Finland; (b) Karlsruhe, Germany; (c) Park Falls, USA; and (d) Lauder, New 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 ${\mathrm{XCO}}_{\mathrm{2}}$ (solid lines) and the daily averaged simulated ${\mathrm{XCO}}_{\mathrm{2}}$ in the observation space (coloured dots). The bottom panel presents the weekly averaged bias of the simulated ${\mathrm{XCO}}_{\mathrm{2}}$ against the TCCON data (coloured dots) and the smoothed bias (solid lines). Blue represents the free run, while red is for the analysis.

[Figure omitted. See PDF]

Global evaluation of the analysis

In this section we first present the characteristics of the ${\mathrm{XCO}}_{\mathrm{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 ${\mathrm{XCO}}_{\mathrm{2}}$ comparing the ${\mathrm{XCO}}_{\mathrm{2}}$ from the analysis against the ${\mathrm{XCO}}_{\mathrm{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 $\mathit{\delta }$ of $-$0.36 $\mathrm{ppm}$ and a mean absolute error $\mathrm{\Delta }$ of 1.08 $\mathrm{ppm}$ (Table ). However, the individual station bias ${\mathit{\delta }}_k$ spans a range from 2.3 $\mathrm{ppm}$ at Ascension Island (Saint Helena, Ascension and Tristan da Cunha) to $-$2.9 $\mathrm{ppm}$ at Białystok (Poland). The station-to-station bias deviation $\mathit{\sigma }$ 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. a. First, it shows that the initial condition of the free run has a positive bias of about 2 $\mathrm{ppm}$ over the tropical region (region between 23${}^{\circ }$ S and 23${}^{\circ }$ 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 $\mathrm{ppm}$. These results are slightly different from those of , where the model bias was found to be more constant in the tropical region when comparing the background ${\mathrm{CO}}_{\mathrm{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 $\mathrm{XCO}\mathrm{2}$ measurements from the TCCON 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-latitudes (region between 23 and 66${}^{\circ }$ N) and reaches almost 4 $\mathrm{ppm}$ at the latitude of Sodankylä (Finland, 67${}^{\circ }$ N) when compared to the TCCON ${\mathrm{XCO}}_{\mathrm{2}}$. This value is the result of the smoothing process as we do not have data for that period (Fig. a). 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. b and c). The negative bias at northern mid-latitudes remains high during the whole year, with an absolute value generally greater than 1 $\mathrm{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 ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$ sink associated with the growing season starts too early in the season .

The precision $\mathit{\pi }$ of the free run measured by the average scatter between the simulation and the TCCON data is 1.4 $\mathrm{ppm}$ (Table ). Similarly to the bias, the scatter varies in time and space as highlighted by the Hovmöller diagram of the scatter (Fig. a). The scatter has its highest values of more than 1 $\mathrm{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 $\mathrm{ppm}$ in 36 h, whereas the measurements show a decrease of about 4 $\mathrm{ppm}$ (Fig. a). 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 ${\mathrm{CO}}_{\mathrm{2}}$.

Analysis vs. free run

To assess the impact of the assimilation of the MACC GOSAT BESD ${\mathrm{XCO}}_{\mathrm{2}}$, we compared the evolution of ${\mathrm{XCO}}_{\mathrm{2}}$ from the analysis with ${\mathrm{XCO}}_{\mathrm{2}}$ from the free simulation. Figure  presents the Hovmöller diagram (time vs. latitude) of this difference. It shows that the first region where the analysis impacts ${\mathrm{XCO}}_{\mathrm{2}}$ is the tropics. There, compared to the free run, the analysis continuously decreases ${\mathrm{XCO}}_{\mathrm{2}}$ by up to 1 $\mathrm{ppm}$ in June and by more than 2 $\mathrm{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 ${\mathrm{XCO}}_{\mathrm{2}}$ over the southern extra tropics (region between 23 and 66${}^{\circ }$ S) when compared to the free run (Fig. ). The decrease extends to the southern high latitudes ($\ge$ 66${}^{\circ }$ S) even when no GOSAT data were assimilated in this region. This decrease results mainly from the transport of ${\mathrm{CO}}_{\mathrm{2}}$ from the equatorial region and southern mid-latitudes towards southern high latitudes. Unfortunately, there are no independent ${\mathrm{XCO}}_{\mathrm{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${}^{\circ }$ N, the analysis has higher values of ${\mathrm{XCO}}_{\mathrm{2}}$ than the free run, with a difference of more than 2 $\mathrm{ppm}$ during the northern summer. Again, the assimilation of the GOSAT data improves the simulated ${\mathrm{XCO}}_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{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${}^{\circ }$ N.

Analysis vs. TCCON data

When compared with the TCCON data, the GOSAT BESD ${\mathrm{XCO}}_{\mathrm{2}}$ analysis has an offset $\mathit{\delta }$ of $-$0.34 $\mathrm{ppm}$ and a mean absolute error $\mathrm{\Delta }$ of 0.57 $\mathrm{ppm}$ (Table ). The offset is similar to that of the free run ($-$0.36 $\mathrm{ppm}$), but the mean absolute error is improved (1.08 $\mathrm{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 $\mathrm{ppm}{\mathrm{yr}}^{-\mathrm{1}}$ for Lauder (New Zealand) (Table S1 of the Supplement), and it improves to 0.47 $\mathrm{ppm}{\mathrm{yr}}^{-\mathrm{1}}$ for the analysis (Table S2 and Fig. c).

By increasing ${\mathrm{XCO}}_{\mathrm{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 $\mathrm{ppm}$ (Fig. b). This could be explained by the fact that we correct the atmospheric state of ${\mathrm{CO}}_{\mathrm{2}}$ and not the ${\mathrm{CO}}_{\mathrm{2}}$ fluxes. During the seasons when the ${\mathrm{CO}}_{\mathrm{2}}$ fluxes are the main driver of the atmospheric ${\mathrm{CO}}_{\mathrm{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 $\mathit{\sigma }$ that is largely reduced compared to the free run with a value of 0.61 $\mathrm{ppm}$ against 1.27 $\mathrm{ppm}$ (Table ). The assimilation of the MACC GOSAT BESD ${\mathrm{XCO}}_{\mathrm{2}}$ thus helps to significantly improve the accuracy of the model. The assimilation also helps improve the precision $\mathit{\pi }$, with the mean scatter improved by 15 $\mathrm{\%}$, reduced to a value of 1.22 $\mathrm{ppm}$. The scatter of the analysis is reduced for all TCCON stations compared to the free run except for Garmisch (Germany), 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. ). In particular, the analysis shows less spurious variability than the free run at Sodankylä (Fig. a).

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 TCCON data using a geolocation criterion of 5${}^{\circ }$ in space and a time window of $\pm \mathrm{2}$ $\mathrm{h}$. Before computing the difference between each GOSAT–TCCON pair, following , we added a correction to the GOSAT-retrieved value in order to account for the use of different a priori ${\mathrm{CO}}_{\mathrm{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 procedure 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 ).

Statistics of the ${\mathrm{XCO}}_{\mathrm{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 (${\mathit{\delta }}_k$, in $\mathrm{ppm}$), scatter (${\mathit{\sigma }}_k$, in $\mathrm{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 $\mathrm{ppm}$) and the mean $r$ (last three rows). The second column ($N$) is the number of data points used for computing the statistics.

		MACC GOSAT data set				Analysis
Site	$N$	Bias	Scatter	$r$		Bias	Scatter	$r$
Sodankylä	90	$-$0.26	4.50	0.39		0.24	1.41	0.92
Białystok	58	$-$0.28	3.45	0.32		1.06	1.99	0.17
Bremen	41	1.19	2.34	0.53		0.54	0.86	0.81
Karlsruhe	91	1.45	2.74	0.52		0.89	0.74	0.88
Orléans	52	0.20	2.44	0.34		1.29	0.57	0.84
Garmisch	76	1.64	3.10	0.55		1.17	1.06	0.77
Park Falls	63	1.50	3.22	0.71		$-$0.08	1.03	0.95
Four Corners	102	$-$0.00	3.79	0.64		0.65	0.81	0.89
Lamont	340	$-$1.01	4.05	0.57		0.05	1.01	0.91
Saga	61	0.40	2.95	0.76		0.14	0.88	0.90
Darwin	234	$-$1.27	3.37	0.42		$-$0.11	0.81	0.84
Wollongong	221	$-$3.03	3.86	0.31		$-$1.54	1.07	0.74
Mean	12	0.04	3.32	0.50		0.36	1.02	0.80
MAE	12	1.02	–	–		0.65	–	–
Deviation	12	1.31	–	–		0.74	–	–

Hovmöller diagram (latitude vs. time) of the difference in ppm (negative/positive in blue/red) between ${\mathrm{XCO}}_{\mathrm{2}}$ from the analysis and from the free-run simulation, from 1 January to 31 December 2013. The horizontal dotted lines represent the latitude of the northernmost and the southernmost TCCON station, respectively. The grey shaded areas are where GOSAT does not provide observations.

[Figure omitted. See PDF]

For each GOSAT–TCCON pair, we extracted the ${\mathrm{CO}}_{\mathrm{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  and ). 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 $\mathrm{\Delta }$ than the one from the MACC GOSAT BESD data (0.65 $\mathrm{ppm}$ vs. 1 $\mathrm{ppm}$, Table ), a station-to-station bias deviation $\mathit{\sigma }$ almost half of the one from GOSAT data (0.7 $\mathrm{ppm}$ vs. 1.3 $\mathrm{ppm}$) and has an improved precision $\mathit{\pi }$ (1 $\mathrm{ppm}$ vs. 3.3 $\mathrm{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 , 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 $\sim 0.4$ $\mathrm{ppm}$ and a precision of $\sim \mathrm{2}$ $\mathrm{ppm}$.

The better precision (lower value of $\mathit{\pi }$) and the lower value of the mean absolute error $\mathrm{\Delta }$ and station-to-station bias deviation $\mathit{\sigma }$ 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 ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{XCO}}_{\mathrm{2}}$ dropped from 398.62 $\mathrm{ppm}$ at 08:15 $\mathrm{LT}$ (local time) to 395.97 $\mathrm{ppm}$ at 12:53 $\mathrm{LT}$ (Fig. , top panel). This sudden decrease of 2.65 $\mathrm{ppm}$ in less than 5 $\mathrm{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 $\mathrm{hPa}$ (Fig. , lower panel).

The free run is able to capture the sudden decrease in ${\mathrm{XCO}}_{\mathrm{2}}$, highlighting the skill of the model for such a situation (Fig. , 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 ${\mathrm{CO}}_{\mathrm{2}}$ relative to the background (Figs. e and f). When it moves towards Park Falls, it results in decreasing ${\mathrm{XCO}}_{\mathrm{2}}$ as observed and simulated, but the decrease in the free run is too strong by 2 to 3 $\mathrm{ppm}$ compared to the measurements.

Situation over Park Falls (USA) between 30 May and 2 June. Top panel: evolution of ${\mathrm{XCO}}_{\mathrm{2}}$ (in $\mathrm{ppm}$) from hourly averaged TCCON data (black dots), the free run (blue line and dots) and the analysis (red line and dots). The dots are the values of the model in the observation space. Lower panel: evolution of the mean sea level pressure (in $\mathrm{hPa}$, black line) and the temperature at 500 $\mathrm{hPa}$ (in $\mathrm{K}$, magenta line). The vertical dotted lines represent 31 May, at 00:00 $\mathrm{UTC}$ and at 12:00 $\mathrm{UTC}$, and 1 June, at 00:00 $\mathrm{UTC}$.

[Figure omitted. See PDF]

Situation around Park Falls (black triangle), Wisconsin, USA, end of May 2013. (a) Average increment in terms of ${\mathrm{XCO}}_{\mathrm{2}}$ (in ppm, negative/positive in blue/red) on 30 May 2013 (contours) and location of the GOSAT measurements during this day (black rectangles). (b, c, d) ${\mathrm{XCO}}_{\mathrm{2}}$ (in ppm) on 31 May at 00:00 $\mathrm{UTC}$, at 12:00 $\mathrm{UTC}$ and on 1 June at 00:00 $\mathrm{UTC}$, respectively, from the analysis. (e, f) ${\mathrm{XCO}}_{\mathrm{2}}$ (in ppm) on 31 May at 12:00 $\mathrm{UTC}$ and on 1 June at 00:00 $\mathrm{UTC}$ from the free run (below/above background value in blue/red). For (b) to (f) the dark contours are the values of the geopotential at 500 $\mathrm{hPa}$.

[Figure omitted. See PDF]

We investigated whether the assimilation of the GOSAT data helps improve the simulated evolution of the ${\mathrm{CO}}_{\mathrm{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. a). These measurements have the effect of increasing the ${\mathrm{XCO}}_{\mathrm{2}}$ in this region (Fig. b–d). The cold air mass is then richer in ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{XCO}}_{\mathrm{2}}$ horizontal gradient.

The ${\mathrm{XCO}}_{\mathrm{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 ${\mathrm{XCO}}_{\mathrm{2}}$ decrease.

Anomaly correlation coefficient (ACC) of the forecast compared to its own analysis as a function of the forecast lead time and for each month: (a) global ACC, (b) ACC for the Northern Hemisphere (20–90${}^{\circ }$ N), (c) ACC for the tropics (20${}^{\circ }$ S–20${}^{\circ }$ N), and (d) ACC for the Southern Hemisphere (90–20${}^{\circ }$ S). Each month is represented by a different colour (see inset legends).

[Figure omitted. See PDF]

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 10-day 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 ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{XCO}}_{\mathrm{2}}$ (see Appendix  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 $\mathrm{\%}$ corresponds to forecasts for which the error is the same as for a forecast based on a climatological average. An ACC of about 80 $\mathrm{\%}$ 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 $\mathrm{h}$ up to 10 days. We found that the ACC is globally more than 90 $\mathrm{\%}$ for day 3 and almost always more than 85 $\mathrm{\%}$ for day 5 for each single month (Fig. a). This means that the forecast for today based on the analysis of 5 days ago shows the same large-scale synoptic ${\mathrm{XCO}}_{\mathrm{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. b to d). This difference between the two hemispheres may reflect the fact that the ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{XCO}}_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$ fluxes are an even more important driver of the ${\mathrm{CO}}_{\mathrm{2}}$ concentration than the initial ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$ in the atmosphere in the model over the Northern Hemisphere (as confirmed by Fig. a). This induces a negative bias of the ${\mathrm{CO}}_{\mathrm{2}}$ surface fluxes in the model. In the meantime, the analysis increases the ${\mathrm{CO}}_{\mathrm{2}}$ concentration helped by the GOSAT BESD data (Fig. ). However, the next 12 h short-term forecast (used as the background for the next analysis) will not increase the ${\mathrm{CO}}_{\mathrm{2}}$ concentration enough due to the negative bias of the ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$ fluxes are important and there is a positive bias in the fluxes (too much release or too little uptake of ${\mathrm{CO}}_{\mathrm{2}}$ in the atmosphere) in the model. For these situations when the ${\mathrm{CO}}_{\mathrm{2}}$ fluxes are the main driver of the atmospheric ${\mathrm{CO}}_{\mathrm{2}}$, the only solution to improve the skill would be to optimize the ${\mathrm{CO}}_{\mathrm{2}}$ fluxes together with the ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$ instead of the surface fluxes in order to constrain the atmospheric ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$, model errors are related to potentially inaccurate surface fluxes or transport.

This article demonstrates the benefit of the assimilation of ${\mathrm{XCO}}_{\mathrm{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 ${\mathrm{XCO}}_{\mathrm{2}}$ provides a ${\mathrm{CO}}_{\mathrm{2}}$ analysis that was compared to a free-run forecast where the ${\mathrm{CO}}_{\mathrm{2}}$ concentration is not constrained by any ${\mathrm{CO}}_{\mathrm{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 and mainly associated with biogenic fluxes.

The analysis significantly reduces these biases without completely removing them, with a remaining mean offset of $-$0.34 $\mathrm{ppm}$ and a mean absolute error of 0.57 $\mathrm{ppm}$ compared to the TCCON data. However, the accuracy estimated with the station-to-station bias deviation is 0.61 $\mathrm{ppm}$. This represents a large improvement compared to the free run, for which the accuracy is 1.27 $\mathrm{ppm}$. The precision of the analysis estimated with the mean scatter is 1.22 $\mathrm{ppm}$, slightly better than for the free run with a value of 1.43 $\mathrm{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 $\mathrm{ppm}$ compared to 1.3 $\mathrm{ppm}$). The precision is much better for the analysis, with a scatter of 1 $\mathrm{ppm}$, while the assimilated data have a scatter of 3.3 $\mathrm{ppm}$. The precision of the analysis is also better than the documented precision of other GOSAT ${\mathrm{XCO}}_{\mathrm{2}}$ products. The precision of the NIES product extracted from is 1.8 $\mathrm{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 $\mathrm{ppm}$ . The ${\mathrm{CO}}_{\mathrm{2}}$ analysis is consequently an alternative to the standard ${\mathrm{XCO}}_{\mathrm{2}}$ GOSAT products as it provides a lower or similar station-to-station bias deviation and a better-precision ${\mathrm{XCO}}_{\mathrm{2}}$ product compared to TCCON. Moreover, it has a uniform spatio-temporal resolution.

The pre-operational CAMS ${\mathrm{CO}}_{\mathrm{2}}$ analysis is similar to the analysis presented in this paper, having nevertheless a higher horizontal resolution (TL511 on a reduced Gaussian grid, $\sim 40\mathrm{km}\times 40\mathrm{km}$), and a higher vertical resolution with 137 vertical levels. It currently assimilates the most recent version of the GOSAT BESD data presented by in near-real time. These data have an improved bias deviation ($\sim$ 0.4 $\mathrm{ppm}$) and an improved precision ($\sim$ 2 $\mathrm{ppm}$) compared to those used in this study. The near-real-time CAMS ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$ model used in this study is its ability to represent ${\mathrm{CO}}_{\mathrm{2}}$ variations associated with synoptic weather systems . By correcting the large-scale ${\mathrm{XCO}}_{\mathrm{2}}$ patterns and removing part of the model bias, we showed with a case study that the analysis is able to better represent the ${\mathrm{CO}}_{\mathrm{2}}$ variations associated with these situations. The variations in the atmospheric reservoir of ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{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\mathrm{km}\times 16\mathrm{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 $\mathrm{\%}$, respectively. This means that we are providing a ${\mathrm{CO}}_{\mathrm{2}}$ forecast with accurate synoptic features for today. With a good representation of the variability and a bias mostly under 1 $\mathrm{ppm}$, the CAMS atmospheric ${\mathrm{CO}}_{\mathrm{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).

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 ${\widehat{c}}^o$, the TCCON-retrieved ${\mathrm{XCO}}_{\mathrm{2}}$ to ${\mathit{x}}^{\mathrm{t}}$, the true (unknown) ${\mathrm{CO}}_{\mathrm{2}}$ profile , $${\widehat{c}}^o=c^b+{\mathit{a}}^{\mathrm{T}}\left({\mathit{x}}^{\mathrm{t}},-,{\mathit{x}}^b\right)+\mathit{\epsilon },$$ where ${\mathit{x}}^b$ is an a priori profile of ${\mathrm{CO}}_{\mathrm{2}}$, $\mathit{a}$ is a vector resulting from the product of the averaging kernel matrix with a dry-pressure weighting function vector (for the vertical integration), $c^b$ is the column-averaged mixing ratio computed from ${\mathit{x}}^b$, and $\mathit{\epsilon }$ 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 $\mathit{x}$ is converted to a column-averaged mixing ratio $\widehat{c}$ by $$\widehat{c}=c^b+{\mathit{a}}^{\mathrm{T}}\left(\mathit{x},-,{\mathit{x}}^b\right).$$

The comparison between the simulation and TCCON occurs in the observation space with the difference between the model column-averaged mixing ratio $\widehat{c}$ of Eq. () and the TCCON column-averaged mixing ratio ${\widehat{c}}^o$ of Eq. (), $$\widehat{c}-{\widehat{c}}^o={\mathit{a}}^{\mathrm{T}}\left(\mathit{x},-,{\mathit{x}}^{\mathrm{t}}\right)-\mathit{\epsilon }.$$

Let us define $\mathit{\eta }={\mathit{a}}^{\mathrm{T}}(\mathit{x}-{\mathit{x}}^{\mathrm{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 ${\mathrm{XCO}}_{\mathrm{2}}$ data for the analysis. The difference between the smooth model column-averaged mixing ratio $\widehat{c}$ and the TCCON column-averaged mixing ratio ${\widehat{c}}^o$ is, therefore, the sum of the model error $\mathit{\eta }$ and the error in the retrieved column-averaged mixing ratio $\mathit{\epsilon }$.

To compute the model column-averaged mixing ratio $\widehat{c}$ of Eq. () 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. ().

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. ). 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. ).

Time smoothing

For each TCCON station $k$ and for each week $w^l$ for $l\in \left[\mathrm{1},,,52\right]$, we compute the mean difference ${\mathit{\delta }}_k^l$ and the standard deviation of the difference ${\mathit{\sigma }}_k^l$ 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 ${\stackrel{\mathrm{̃}}{b}}_k\left(t\right)$ that combines a linear growth and a harmonic component, $$\begin{array}{cc}& {\stackrel{\mathrm{̃}}{b}}_k\left(t\right)=a_{k}t+b_k+{\mathit{\alpha }}_k\sin \left({\displaystyle \frac{t}{{\mathit{\tau }}_{\mathrm{1}}}},+,{\mathit{\phi }}_k\right)\\ & & +{\mathit{\beta }}_k\sin \left({\displaystyle \frac{t}{{\mathit{\tau }}_{\mathrm{2}}}},+,{\mathit{\phi }}_k\right).\end{array}$$ $a_k$, $b_k$, ${\mathit{\alpha }}_k$, ${\mathit{\beta }}_k$ and ${\mathit{\phi }}_k$ are the parameters of the function ${\stackrel{\mathrm{̃}}{b}}_k\left(t\right)$ obtained by an optimization procedure that minimizes the distance between ${\stackrel{\mathrm{̃}}{b}}_k\left(t\right)$ and the series of ${\mathit{\delta }}_k^l$ for $l\in \left[\mathrm{1},,,52\right]$. ${\mathit{\tau }}_{\mathrm{1}}$ is chosen to be 6 months and ${\mathit{\tau }}_{\mathrm{2}}$ 3 months. The form of the function of Eq. () thus gives a linear growing bias and allows seasonal variations. A similar function is used for the standard deviation.

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 , we then compute for each week $w^l$ the best fit of the interpolated biases with a quadratic function of latitude ${\widehat{b}}^l$, $${\widehat{b}}^l\left(\mathit{\varphi }\right)=a^l{\mathit{\varphi }}^{\mathrm{2}}+b^l\mathit{\varphi }+c^l,$$ where $\mathit{\varphi }$ is the sine of the latitude. $a^l$, $b^l$ and $c^l$ are obtained by an optimization procedure that minimizes the distance between ${\widehat{b}}^l$ and the weekly interpolated biases ${\mathit{\delta }}_k^l$ for $k\in \left[\mathrm{1},,,N\right]$. A similar function is used for the standard deviation.

Discussion

For some stations, the availability of the weekly differences is not uniform in time and the time smoothing of Eq. () provides spurious values. We solved this issue by fixing the coefficient ${\mathit{\alpha }}_k$ to a zero value (See Table S1).

With a root mean square error (RMSE) mostly under 0.7 $\mathrm{ppm}$ and a correlation mostly over 0.8, the smoothed bias matches well with the weekly bias (Table S1). The Hovmöller diagram (Fig. ) 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).

Anomaly correlation coefficient

The anomaly correlation coefficient (ACC) between the forecast $f$ and the analysis $a$ is computed using the climatology $c$ by $$\text{ACC}={\displaystyle \frac{\stackrel{\mathrm{‾}}{(f-c)(a-c)}}{\sqrt{\stackrel{\mathrm{‾}}{(f-c{)}^{\mathrm{2}}}\stackrel{\mathrm{‾}}{(a-c{)}^{\mathrm{2}}}}}},$$ where the overline is the spatial and temporal average. For example, for the forecast range 24 h, we take the ${\mathrm{XCO}}_{\mathrm{2}}$ fields from all the 24 h forecasts for a given month, all the analyses valid for the same time, and a fixed climatology for this month.

The climatology is based on a free-run simulation using the optimized ${\mathrm{CO}}_{\mathrm{2}}$ surface fluxes from , which simulated the years from 2003 to 2012. For each month, we compute the average over the 10 years of the simulation, rescaling the mean so that the mean is the same as for the analysis, avoiding with this procedure the issue of the increase in ${\mathrm{CO}}_{\mathrm{2}}$ over time. The two-dimensional climatology field for ${\mathrm{XCO}}_{\mathrm{2}}$ for the month $m$ is $$\begin{array}{cc}& c\left(m\right)={\displaystyle \frac{\mathrm{1}}{10}}\sum _{y=2003}^{2012}{\displaystyle \frac{\mathrm{1}}{n(y,m)}}\sum _{d=\mathrm{1}}^{n(y,m)}\\ & & \left[\mathrm{\Sigma },\left(y,,,m,,,d\right),-,\stackrel{\mathrm{‾}}{\mathrm{\Sigma }},\left(y,,,m,,,d\right)\right]+{\stackrel{\mathrm{‾}}{\mathrm{\Sigma }}}_{\text{an}}\left(m\right),\end{array}$$ where $y$ is the year, $n$ the number of days for the year $y$ and the month $m$, $d$ is an index for the day, $\mathrm{\Sigma }\left(y,,,m,,,d\right)$ is the ${\mathrm{XCO}}_{\mathrm{2}}$ field from the simulation for the year $y$, the month $m$ and the day $d$, and $\stackrel{\mathrm{‾}}{\mathrm{\Sigma }}$ is a spatial average of $\mathrm{\Sigma }$ and ${\stackrel{\mathrm{‾}}{\mathrm{\Sigma }}}_{\text{an}}\left(m\right)$ is the spatial and temporal average of the ${\mathrm{XCO}}_{\mathrm{2}}$ fields from the analysis for the month $m$ (and the year 2013).

The Supplement related to this article is available online at doi:10.5194/acp-16-1653-2016-supplement.

S. Massart designed and carried out the experiments with the help of A. Agustí-Panareda and advice from F. Chevallier, J. Heymann and M. Buchwitz. J. Heymann, M. Reuter, M. Hilker, M. Buchwitz and J. P. Burrows were responsible for the design and operation of the BESD GOSAT ${\mathrm{XCO}}_{\mathrm{2}}$ retrieval algorithm. S. Massart prepared the manuscript with contributions from A. Agustí-Panareda, J. Heymann, M. Buchwitz, F. Chevallier, M. Reuter, M. Hilker, J.P. Burrows, D. G. Feist, and F. Hase. N. M. Deutscher and R. Sussmann contributed to the ACP version of the paper. F. Desmet is the co-investigator of the La Réunion TCCON station. N. M. Deutscher and C. Petri are responsible for the Białystok, Bremen and Orléans TCCON data. M. Dubey is the PI of Four Corners TCCON station. D. G. Feist is the PI of the Ascension TCCON station. D. W. T. Griffith and V. Velazco are the PIs of Darwin and Wollongong stations. F. Hase is the PI of the Karlsruhe TCCON station. R. Kivi is the PI of the Sodankylä TCCON station. M. Schneider is the PI of Izaña TCCON station. R. Sussmann is the PI of the Garmisch TCCON station.

Acknowledgements

This study was funded by the European Commission under the European Union's Horizon 2020 programme. The development of the GOSAT BESD algorithm received funding from the European Space Agency (ESA) Greenhouse Gases Climate Change Initiative (GHG-CCI). TCCON data were obtained from the TCCON Data Archive, hosted by the Carbon Dioxide Information Analysis Center (CDIAC) – http://tccon.ornl.gov/. Garmisch work was funded in part via the ESA GHG-CCI project. Four Corners TCCON was funded by LANL's LDRD programme. Darwin and Wollongong TCCON measurements are funded by NASA grants NAG5-12247 and NNG05-GD07G and the Australian Research Council grants DP140101552, DP110103118, DP0879468, LE0668470 and LP0562346. We are grateful to the DOE ARM programme for technical support in Darwin, and Clare Murphy, Nicholas Jones and others for support in Wollongong. TCCON measurements in Białystok and Orléans are supported by ICOS-INWIRE, InGOS and the Senate of Bremen. N. Deutscher is supported by an ARC-DECRA fellowship, DE140100178. The authors are grateful to Marijana Crepulja for the acquisition of the BESD GOSAT data at ECMWF and the preparation of the data for the assimilation. The authors would like to acknowledge Paul Wennberg, PI of the Lamont and Park Falls TCCON stations. Finally, we would like to express our great appreciation to William Lahoz, editor of this paper, for his useful comments during the revision process. Edited by: W. Lahoz