Forecasting global atmospheric CO 2

Abstract. A new global atmospheric carbon dioxide (CO2) real-time forecast is now available as part of the pre-operational Monitoring of Atmospheric Composition and Climate – Interim Implementation (MACC-II) service using the infrastructure of the European Centre for Medium-Range Weather Forecasts (ECMWF) Integrated Forecasting System (IFS). One of the strengths of the CO2 forecasting system is that the land surface, including vegetation CO2 fluxes, is modelled online within the IFS. Other CO2 fluxes are prescribed from inventories and from off-line statistical and physical models. The CO2 forecast also benefits from the transport modelling from a state-of-the-art numerical weather prediction (NWP) system initialized daily with a wealth of meteorological observations. This paper describes the capability of the forecast in modelling the variability of CO2 on different temporal and spatial scales compared to observations. The modulation of the amplitude of the CO2 diurnal cycle by near-surface winds and boundary layer height is generally well represented in the forecast. The CO2 forecast also has high skill in simulating day-to-day synoptic variability. In the atmospheric boundary layer, this skill is significantly enhanced by modelling the day-to-day variability of the CO2 fluxes from vegetation compared to using equivalent monthly mean fluxes with a diurnal cycle. However, biases in the modelled CO2 fluxes also lead to accumulating errors in the CO2 forecast. These biases vary with season with an underestimation of the amplitude of the seasonal cycle both for the CO2 fluxes compared to total optimized fluxes and the atmospheric CO2 compared to observations. The largest biases in the atmospheric CO2 forecast are found in spring, corresponding to the onset of the growing season in the Northern Hemisphere. In the future, the forecast will be re-initialized regularly with atmospheric CO2 analyses based on the assimilation of CO2 products retrieved from satellite measurements and CO2 in situ observations, as they become available in near-real time. In this way, the accumulation of errors in the atmospheric CO2 forecast will be reduced. Improvements in the CO2 forecast are also expected with the continuous developments in the operational IFS.


Introduction
Atmospheric composition monitoring was integrated in the numerical weather prediction framework (NWP) at the European Centre for Medium-Range Weather Forecasts (ECMWF) as part of the Global and regional Earth-System Monitoring using Satellite and in situ data (GEMS) and the Monitoring of Atmospheric Composition and Climate (MACC) projects (Hollingsworth et al., 2008). The resulting global forecasting system of atmospheric composition benefits from the existing operational infra-structure for weather forecasting, satellite data assimilation and high performance computing at ECMWF. Until recently, only forecasts of reactive gases and aerosols were provided in near-real time on a routine basis Morcrette et al., 2009) Published by Copernicus Publications on behalf of the European Geosciences Union.

11960
A. Agustí-Panareda et al.: Global CO 2 forecast as part of the Copernicus European programme, formerly called GMES (Global Monitoring for Environment and Security). The reasons for not having carbon dioxide (CO 2 ) stemmed from the challenges associated with modelling the CO 2 fluxes and the relatively small signals characterizing CO 2 variability making the accuracy requirements for the model simulations more stringent than for other trace gases. The recent addition of the CTESSEL carbon module in the operational Integrated Forecasting System (IFS) at ECMWF (Boussetta et al., 2013a) has now also made feasible the delivery of atmospheric CO 2 forecasts in real time. Although the forecast is currently not initialized with a CO 2 analysis because of the lack of CO 2 observations with global coverage in near-real time, it relies heavily on a wealth of meteorological observations for initializing the meteorology and transport. Moreover, we expect that in the near future there will be satellite retrievals of CO 2 from the GOSAT (Greenhouse Gases Observing Satellite)(www.gosat.nies.go.jp) and the OCO-2 (Orbiting Carbon Observatory) (http://oco.jpl. nasa.gov) available a few days behind real time. These CO 2 satellite products will be assimilated to produce CO 2 analyses also in near-real time. It is worth noting that the CO 2 retrievals provide averaged column information of CO 2 and only for sunlit clear-sky conditions. Therefore, they cannot provide information on the CO 2 vertical distribution, neither at nighttime and during winter time at high latitudes, nor on the CO 2 anomalies associated with cloudy regions within convective and synoptic weather systems. Thus, the CO 2 forecast model will be crucial in filling this information gap during the data assimilation process. Indeed, the main use of the forecast is to support the data assimilation of CO 2 observations. Because the data assimilation window used in the IFS is 12 h, the main requirement for the CO 2 forecast is to have skill in the simulation of the CO 2 variability on short timescales, e.g. diurnal and synoptic scales. The errors in the forecast will influence the quality of the resulting CO 2 analysis. For this reason, the evaluation of the CO 2 forecast errors is also very important for the analysis. The in situ observations at the surface are very valuable not only for evaluation purposes, but they have the potential to provide complementary information to the CO 2 satellite products for the CO 2 analysis. The continuous in situ observations are much more accurate than the satellite data, therefore providing a reference for correcting biases close to the surface. Although they have a sparser spatial coverage than satellite measurements, they have a much better temporal coverage at high latitudes, during cloudy conditions and at nighttime.
The atmospheric CO 2 variability results mainly from a strong synergy between surface fluxes and atmospheric transport. The advection of CO 2 across meridional gradients associated with large-scale flux patterns dominates the variability in the free troposphere, whereas local fluxes also play a role in the variability of atmospheric CO 2 close to the surface, i.e. within the atmospheric boundary layer (Keppel-Aleks et al., 2011, 2012. Modelling the spatial and temporal CO 2 variability is a challenging task. The difficulties arise from uncertainties in the modelling of both the sources/sinks (le Quéré et al., 2009) and transport Patra et al., 2008).
Globally, the CO 2 variability on timescales ranging from diurnal, seasonal, to interannual is dominated by the terrestrial biogenic fluxes (Geels et al., 2004). The challenge of modelling the terrestrial biogenic fluxes comes from high spatial heterogeneity of the land surface and complex processes with large uncertainties. Some of these uncertainties stem from a lack of observational data with sufficient global coverage to characterize all the variability in space and time associated with vegetation and carbon pools. At the same time, the biospheric fluxes are strongly influenced by climate variability (Keeling et al., 1995). Therefore, the timely availability of accurate meteorological data sets is also crucial. The recent development of the CTESSEL carbon module within the IFS takes advantage of accurate real-time climate forcing in order to provide online terrestrial biogenic fluxes also in real time.
The online computation of terrestrial biogenic fluxes and transport -both forced and initialized by NWP analysesis key to ensure consistency in the coupling between fluxes and transport. An example of the importance of this consistency is the passage of mid-latitude frontal cyclones. The change in radiation associated with the frontal cloud reduces the photosynthetic CO 2 uptake which results in a substantial increase in atmospheric CO 2 (∼ 10 ppm) near the surface, as respiration continues to emit CO 2 (Chan et al., 2004). This high-CO 2 anomaly can then be transported by frontal ascent to the mid and upper troposphere. This coupling between fluxes and transport also works on a seasonal scale. Meridional transport by mid-latitude frontal cyclones reduces/amplifies the seasonal cycle at mid/high latitudes (Parazoo et al., 2011). On diurnal scales and seasonal timescales, there is a covariance between turbulent mixing in the planetary boundary layer and terrestrial biogenic fluxes known as the rectifier effect (Denning et al., 1999).
In addition to the modelling challenges, the availability of CO 2 observations is central to be able to provide optimal estimates of CO 2 concentrations and fluxes, as well as error estimates of the CO 2 model forecasts. So far, the most accurate CO 2 observations are from in situ measurements close to the surface. In the past, these have been available with a delay of 1 to 2 years. This long delay in the availability of observations -combined with the large uncertainties in modelling of fluxes, their forcings, and the transport model -has hindered the task of providing CO 2 information in a timely manner. However, the Integrated Carbon Observation System (ICOS) observing network recently started to provide continuous in situ CO 2 observations with a 1-day lag as part of their preoperational phase. Currently, there are seven stations in the pre-operational network. Some of these stations are sampling baseline air, and therefore allow a continuous monitoring of the background bias in the CO 2 forecast.
Atmos. Chem. Phys., 14, 11959-11983, 2014 www.atmos-chem-phys.net/14/11959/2014/ The aim of this paper is to document the capabilities and limitations of this real-time CO 2 forecast, currently available with a 5 day lead time. This is done by comparing CO 2 hindcasts -i.e model simulations for the past 10 years using the same configuration as the real-time CO 2 forecast -with a wide range of independent observations, thus giving an assessment of the representation of the CO 2 spatial and temporal variability at different scales. Furthermore, the continuous automated monitoring of the atmospheric CO 2 forecast with ICOS observations is also shown. This evaluation supports the ongoing monitoring of the model errors. It is also the first step towards being able to assimilate CO 2 observations in near-real time.
The paper is structured as follows. The description of the model CO 2 fluxes and transport is presented in Sect. 2. The evaluation of the CO 2 hindcasts is done in Sect. 3 by using observations from the Integrated Carbon Observation System (ICOS), the Total Carbon Column Observing Network (TCCON), the National Oceanic and Atmospheric Administration (NOAA) networks and the HIAPER Pole to Pole Observations (HIPPO) field experiment. The CO 2 hindcast performance is discussed in Sect. 4, highlighting future work to reduce the errors as part the operational upgrades of the system. Finally, Sect. 5 recaps on the CO 2 forecast capabilities and possible applications.

Forecast configuration and model description
This section presents the CO 2 forecast set-up, including a description of the transport and flux components in the model. The CO 2 modelling is done within the NWP framework, using the IFS model from ECMWF. Both transport and terrestrial biogenic carbon fluxes are computed online and other prescribed fluxes are read from inventories. This ensures a consistency between flux resolution and transport resolution and it also allows a full coupling between meteorological forcing of biogenic fluxes and transport. A description of the main features of the IFS transport are provided in Sect. 2.1. Section 2.2 describes the different fluxes included in the model in more detail.
In order to be able to evaluate the CO 2 forecast over different time scales, yearly CO 2 hindcasts were performed from 2003 to 2012. The hindcasts are made of 24 h forecasts and the meteorological fields are initialized at the beginning of each forecast with ECMWF operational analyses (Rabier et al., 2000;Janisková and Lopez, 2013). Atmospheric CO 2 is initialized on 1 January each year, using the dry molar fraction fields from the optimized fluxes provided by the MACC flux inversion system . In the subsequent forecasts, the atmospheric CO 2 is cycled from one 24 h forecast to the next one, being free to evolve in the model without constraints from CO 2 observations.
In this paper, we present results from the hindcasts with a horizontal resolution corresponding to approximately 80 km and 60 vertical levels, which is the same resolution as the current ECMWF re-analysis (ERA-Interim). This resolution is at the higher end of commonly used resolutions in global chemical transport models (CTMs) .

Transport
The modelling of the transport is performed by the IFS model operational at ECMWF. The model advection is computed by a semi-Lagrangian scheme (Hortal, 2002;Untch and Hortal, 2006). Because it is not mass conserving by default, a proportional global mass fixer is used to ensure the total global budget in the model is conserved from one model time step to the next during advection. The global proportional mass fixer consists of re-scaling the 3-D field of the atmospheric CO 2 mixing ratio by using a global scaling factor. This factor is obtained by dividing the globally integrated atmospheric CO 2 mass before the semi-Lagrangian advection in the model by the one after the advection. The boundary layer mixing is described in Beljaars and Viterbo (1998) and Koehler et al. (2011). The convection scheme is based on Tiedtke (1989) (see Bechtold et al., 2008, for further details). Full documentation on the IFS can be found in www.ecmwf.int/research/ifsdocs. Note that the system presented in this paper is based on model version CY38R1, which was operational from 19 June 2012 to 25 June 2013.
Results from a recent TRANSCOM model intercomparison experiment show that the IFS has relatively accurate representation of the large-scale/inter-hemispheric transport, vertical profiles  and convective uplift , with comparable skill to other CTMs participating in the TRANSCOM study, e.g. GEOSChem, PCTM and TM5. The CO 2 and SF 6 diurnal amplitudes which are largely controlled by the boundary layer mixing were also assessed by Law et al. (2008). Their study found that the IFS was one of the models that simulated the diurnal cycles closer to those observed. Higher horizontal resolution with respect to other CTMs was found to be a contributing factor.
It is worth noting that the NWP analysis of meteorological fields is one of the main elements determining the quality of the transport. Locatelli et al. (2013) found that methane time series simulated by IFS using ECMWF meteorological reanalysis were highly correlated to those simulated by TM5 also using the same re-analysis; whereas the average correlation of IFS with other models using different meteorological analysis was lower.
Finally, the IFS provides one of the best weather forecasts in the medium-range (up to 10-days lead time) based on NWP model intercomparison of skill scores (Richardson et al., 2013). Because the IFS is a world leading state-of-theart NWP model, it is also used as a reference for the development of some CTMs, e.g. TM5 (see Krol et al., 2005

CO 2 fluxes
The CO 2 net ecosystem exchange (NEE) fluxes are from the carbon module of the land surface model in the IFS (CTESSEL) developed as part of the Geoland project (www. geoland2.eu). Because the NEE fluxes are computed online, they are available at the same spatial and temporal resolution as the transport model (∼ 80 km resolution, every 45 min). CTESSEL is a photosynthesis-conductance (A-gs) model based on Calvet et al. (1998Calvet et al. ( , 2004; Calvet (2000) and developed originally by Jacobs et al. (1996). It provides CO 2 fluxes as well as evapotranspiration. However, the evapotranspiration in the IFS is currently still based on the Jarvis approach (Jarvis, 1976) instead of the plant physiological approach of CTESSEL. Despite not having a full coupling between evapotranspiration and CO 2 fluxes, there is some consistency between the two fluxes because they both rely on the same underlying representation of vegetation.
The NEE results from the gross primary production (GPP) and the ecosystem respiration (R eco ) fluxes which are computed independently in the model. The GPP represents the photosynthetic fluxes which are driven by radiation, soil moisture, soil temperature and a prescribed satellite MODIS leaf area index (LAI) fixed monthly climatology (http:// landval.gsfc.nasa.gov/) based on a 9-year averaging process (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008) as described in Boussetta et al. (2013a). The ecosystem respiration is given by empirical formulas driven by soil moisture, soil temperature and snow cover. The model parameters affecting the sensitivity of GPP and R eco to temperature, soil moisture and radiation are listed in Table 2 of Boussetta et al. (2013a).
The meteorological forcing of the fluxes is from the NWP forecast, providing full consistency between variability of the fluxes, the meteorology and the transport processes. Because vegetation growth is represented by an LAI climatology, land use change cannot be represented. There is also no direct representation of the different carbon pools, but a reference respiration parameter for each vegetation type is used to simulate the heterotrophic respiration. The reference value is obtained by optimization with respect to flux measurements for the different vegetation types (see Table 1 in Boussetta et al., 2013a).
There are nine low-vegetation types and six highvegetation types based on the Biosphere-atmosphere transfer scheme (BATS) classification (Dickinson et al., 1986). The NEE flux is an area-fraction weighted sum of the NEE for the dominant high and the dominant low vegetation classes at each grid point. The evaluation of CTESSEL NEE fluxes with observations based on 10 day averaged CO 2 fluxes at 34 sites shows that there is an average correlation of 0.65, and an average bias and root mean square error of −0.1 and 1.7 g C m −2 d −1 , respectively. A more detailed description and evaluation of the CTESSEL GPP, R eco and the resulting NEE fluxes can be found in Boussetta et al. (2013a).
The fire emission flux is from GFAS v1.0 (Kaiser et al. (2012), www.copernicus-atmosphere.eu/about/project_ structure/input_data/d_fire/) which is available one day behind real time. It has a daily temporal resolution and a horizontal resolution of 0.5 • × 0.5 • . The fire fluxes are kept constant throughout the 5-day forecast. The ocean flux is from the Takahashi et al. (2009) climatology with monthly mean fluxes at 4 • × 5 • resolution. The anthropogenic fluxes are annual mean fluxes based on the last available year (2008) of the EDGAR version 4.2 inventory (http://edgar.jrc.ec.europa. eu). In order to account for the increase in the emissions since 2008, the growth in anthropogenic emissions beyond 2008 has been represented using a global re-scaling factor. This is based on estimated anthropogenic CO 2 emission trends of −1.4 and +5.9 % for 2009 and 2010, respectively (Global Carbon Project, www.globalcarbonproject.org), and a climatological trend of +3.1 % for 2011 and 2012. Note that the same climatological trend will be used to extrapolate the anthropogenic fluxes to the present in the operational CO 2 forecast.

Evaluation of CO 2 forecasts
The hindcasts have been evaluated for different periods to assess the global annual budget and its interannual variability  Table 1 show the stations used from each network and their location. HIPPO flight data (http://hippo.ornl.gov/ dataaccess, Wofsy et al., 2012) has also been used to evaluate CO 2 in the free troposphere (Sect. 3.5, see flight tracks in Fig. 1). Vertical profiles from the NOAA Global Monitoring Division (GMD) Carbon Cycle Vertical Profile Network (Tans et al., 1996) have been used to assess the vertical gradients in the model from the lower to the mid troposphere. The computations involved in the processing of the CO 2 hindcast for the comparison with observations are described in the Appendix.

Global CO 2 budget and its interannual variability
The model atmospheric CO 2 growth is the result of the addition of all the fluxes shown in Fig. 2a. The CO 2 fluxes in the model are currently not constrained by atmospheric CO 2 observations. Thus, the budget of the total CO 2 emissionsaffected by all the errors in the CO 2 fluxes -does not match the observed atmospheric growth. This leads to an annual global bias in the modelled atmospheric CO 2 . In the case of optimized fluxes , there is a reasonably good fit between their budget and the observed global growth. Hence, they can be used as a reference, representing a current best estimate for the fluxes at the global scale. Note that the optimized fluxes are not available in near real time because they rely on the highly accurate atmospheric CO 2 flask observations which are currently only provided several months after the date. The annual bias of the model varies from year to year because there are two compensating errors opposing each other. Namely, the underestimation of the NEE sink in the Northern Hemisphere (NH) summer and the underestimation of NEE release in the NH winter by 1 to 2 Gt C month −1 (Fig. 2b) compared to the optimized fluxes of Chevallier et al. (2011). Therefore, the sign of the resulting annual global bias depends on which of these errors dominates when integrated over the year. For instance, in 2010 and 2011 the underestimation of the NEE source is larger than the underestimation of the sink, resulting in a negative global annual bias. Whereas in 2012 the opposite occurs, the underestimation of the sink is larger than that of the source, thus the positive annual global bias. The interannual variability of atmospheric growth is modulated by the NEE interannual variability.
The correlation between the modelled and observed global annual atmospheric growth is 0.74. Although the main contributor to the annual NEE global sink is the NH, the tropics are responsible for its large interannual variability (Fig. 2c). The large error associated with this interannual variability stems from several factors. Namely, the high sensitivity of the biogenic fluxes to climate forcing in the model, combined with large uncertainty in the model parameters, as well as missing and simplified processes in CTESSEL. Moreover, the large gaps in the meteorological observing network in the tropics result in higher errors associated with the climate forcing of the NEE fluxes. Assimilation of satellite products (e.g. soil moisture, LAI and CO 2 ) might help in the evaluation and reduction of these uncertainties and associated errors.
The strong seasonal cycle in the global atmospheric growth (see grey curve in Fig. 2b, defined as the sum of all the surface flux components) comes mainly from the NH mid-latitudes between 30 and 66 • N (Fig. 2d). This suggests that the large underestimation of the global seasonal cycle amplitude is likely associated with errors in midlatitude NEE fluxes. The errors associated with the modelling of the seasonal cycle are examined further in the next section.

CO 2 seasonal cycle
The phase and amplitude of the seasonal cycle of CO 2 are very dependent on latitude. Thus, the model is first evaluated using the NOAA GLOBALVIEW-CO2 (2011) data set which displays the integrated effects of surface CO 2 fluxes over large regions at different latitudinal bands (Fig. 3). At first glance, the annual cycle phase and amplitude and latitude dependency appears to be reasonably represented in the hindcast. However, there are clear discrepancies between the hindcast and GLOBALVIEW-CO2 (2011)  The difference between the GLOBALVIEW product based on observations and model is shown in the right panel. The CO 2 hindcast has been sampled at the same locations as the GLOBALVIEW observations and the same data processing described in Masarie and Tans (1995) has been applied.
cember). Secondly, the onset of the CO 2 sink associated with the growing season starts too early in the hindcast. The sharp CO 2 decrease in mid-latitudes depicted by GLOBALVIEW-CO2 (2011) in June starts in May in the hindcast. This also leads to a longer growing season. The combination of these two factors is consistent with the negative global bias shown in Fig. 2. The GLOBALVIEW-CO2 (2011)  background air from the NOAA/ESRL network, total column measurements from the TCCON network and continuous measurements from the ICOS network. Figure 1 shows the location of the observing stations.
The monthly biases at three continuous ESRL/NOAA background sites (Thoning et al., 2012) confirm that the largest biases are in the NH, as shown by the −10 to −5 ppm bias in the summer months at Barrow, Alaska (Fig. 4a). The negative bias increases in the NH growing season from March to June. This is shown by the differential monthly bias, which depicts how the bias changes with respect to the previous month. The stations in the tropics and South Pole also display mainly negative monthly biases in the background air, with smaller magnitudes, typically between −1 and −2 ppm. Every year, the hindcast is re-initialized with fields from optimized flux simulations constrained with CO 2 observations that convey the atmospheric growth. The differential bias in January each year, thus, depicts the adjustment applied in order to correct for the annual global mean bias in the previous year (see blue dots in Fig. 4). The annual bias in the tropics and South Pole sites is consistent with the bias of the global budget shown in Fig. 2a. The largest interannual variability in the annual bias is also found for the tropical sites. This variability is consistent with that of the bias in the annual global budget. The anomalous 2005 positive annual biases of 1.5 and 2 ppm at the tropical and South Pole sites,respectively, are in line with the 2.5 Gt C annual global bias (equivalent to 1.2 ppm).
The results from the total column evaluation (Fig. 5) are consistent with the findings from the surface measurements and GLOBALVIEW comparisons (Figs. 4 and 3). In Sodankylä (Finland) and Bialystok (Poland) we observe the same underestimation of column CO 2 during NH winter. The hindcast also brings forward the onset of the CO 2 drawdown associated with the growing season by a month. Namely, the observed steep total column CO 2 decrease in June at Sodankylä starts in early May in the model. Similarly, at Bialystok the beginning of the observed CO 2 drawdown is May, whereas the modelled total column CO 2 starts decreasing in April. At Park Falls (Wisconsin, USA) total column CO 2 is underestimated before and after the summertime CO 2 drawdown, and at Lamont (Oklahoma, USA) the CO 2 is only underestimated in winter (January, November and December).
The evaluation of the seasonal cycle based on the ICOS stations ( Fig. 6) is similarly in agreement with previous findings. Ivittut (Greenland) and Puijo (Finland) confirm the underestimation of the winter CO 2 respiration, and the negative bias in Mace Head (Ireland) is also consistent with an underestimation of CO 2 which starts in winter and becomes more pronounced in spring. The CO 2 spikes in Mace Head are associated with specific events influenced by local and nearby continental sources/sinks (Biraud et al., 2002). The background stations of Ivittut and Mace Head have a negative annual bias of ∼ −3 ppm whereas Puijo which is affected by local vegetation fluxes has an annual bias of ∼ −5 ppm. Finally, Lamto (Ivory coast) shows a large positive bias during the dry season when the site is influenced by continental biogenic fluxes, and a small bias during the wet season when the monsoon winds advect background CO 2 from the ocean.

CO 2 synoptic variability
An evaluation of the variability associated with synoptic events is performed at three tall tower sites of the NOAA/ESRL network in continental North America (Ar-  Figure 4. Monthly bias (hindcast -observation) of CO 2 dry molar fraction [ppm] at NOAA/ESRL continuous surface sites sampling background air (green triangles) and differential monthly biases (i.e. difference of monthly bias with respect to previous month) as red triangles from 2003 to 2011. The blue dots highlight the adjustment in CO 2 at the beginning of each year when the model is re-initialized with a simulation from optimized fluxes which has a bias close to zero (see text for details). gyle, Park Falls and West Branch, see Andrews et al., 2014). These sites are directly influenced by local land biospheric fluxes, atmospheric transport and their interaction. The skill in representing the day-to-day variability is assessed for different months in Sect. 3.3.1 and the importance of modelling NEE for the synoptic skill is assessed in Sect. 3.3.2.

Forecast skill of day-to-day CO 2 variability
The synoptic variability is evaluated first by computing the correlation between daily mean atmospheric CO 2 from observations and hindcasts (Table 2) at different sampling levels (Table 3). The correlation coefficients are predominantly higher than 0.5 in the winter months -January, February, November and December -and most sites have values between 0.65 and 0.95. Most of the variability is linked to low pressure systems advecting CO 2 across the large-scale meridional gradient, with a small modulation associated with biogenic fluxes indicated by the very low correlations between atmospheric CO 2 and the modelled NEE fluxes (not shown). In general, the CO 2 hindcast is able to accurately represent the variability associated with the advection by synoptic weather systems.
The spring months -from March to May -display very low or not significant correlations. The large errors in spring (both poor correlations and large biases) are likely associated with modelling errors in the GPP and R eco . Spring is a challenging period for carbon models to model NEE because it is characterized by the transition from predominant respiration in winter to predominant photosynthetic uptake. The timing of this shift in the sign of the daily mean NEE has been analyzed in the model at the two sites where the correlation coefficients are lowest (Park Falls and West Branch). In the model the transition occurs at the beginning of March, which is consistent with the concurrent underestimation of the modelled Atmos. Chem. Phys., 14, 11959-11983, 2014 www.atmos-chem-phys.net/14/11959/2014/   In the summer months, correlation coefficients are mostly above 0.5, with slightly lower values at Argyle, Maine. During summer, the local fluxes and local transport (e.g. the height of the nocturnal boundary layer) have a large influence on the synoptic variability, which is reflected by the higher correlations between atmospheric CO 2 and those parameters (not shown). Local circulations, nocturnal stable boundary layers and the high vegetation activity in the summer are all associated with high uncertainties in the model. The correlations are lower due to the combined effect of the large uncertainties in these local influences.
In autumn, the correlations are higher than in summer. From September to November, both synoptic transport by mid-latitude low pressure systems and biogenic fluxes are important. Moreover, the coupling between the transport and the fluxes is crucial. This is illustrated in Fig. 7 showing the day-to-day variability in tower in situ data from Park Falls in September 2010. The model is able to simulate the peaks of CO 2 on the 7, 11, 21, 23-24 and 29 September 2010, all of them associated with the passage of low pressure systems. The correlation coefficient between observed and modelled CO 2 is 0.81. The modelled and observed CO 2 are similarly correlated with surface pressure (correlation coefficient r = −0.52 and r = −0.56, respectively) and NEE (r = 0.58 and r = 0.57, respectively).
The persistence effect is the main hypothesis to explain the difference in the atmospheric CO 2 errors between spring and autumn. The seasonal cycle amplitude of the NEE budget in CTESSEL is too weak (see Fig. 2b), i.e. respiration/photosynthesis are too weak in the winter/summer. This persistence effect will lead to an early drawdown in spring (due to the winter negative bias), but in autumn the positive bias associated with the weak sink will be compensated by the previous spring negative bias.

Impact of NEE day-to-day variability on the atmospheric CO 2 synoptic forecast skill
The relative importance of the synoptic variability of NEE vs. transport can be assessed by comparing the standard hindcast with a simulation using 3 hourly monthly mean NEE from CTESSEL (i.e. without day-to-day variability) instead of real-time NEE. In order to demonstrate this, it is important to first find observing sites which are systematically affected by both NEE and synoptic advection, and properly represented in the model. The observing station at Park Falls experiences the ideal conditions in September. Both local NEE fluxes and synoptic advection are important for the simulation of the variability of the atmospheric CO 2 there. In addition, the site exhibits a good correlation between the simulated and the observed CO 2 . Figure 8 shows the day-to-day variability of daily minimum, maximum and mean CO 2 at 30 and 396 m level above the surface from the tall tower at Park Falls for the two simulations and observations in September 2010. The observed CO 2 variability is characterized by a trend associated with the seasonal cycle and day-to-day synoptic variability. The variability of the minimum CO 2 during day time is dominated by the trend. Whereas at night time, the CO 2 maximum is modulated by synoptic variations. As expected, the CO 2 day-time trend is present in the hindcast with real-time NEE, but absent in the simulation with the monthly mean NEE. The underestimation of the trend in the hindcast with real-time NEE is consistent with the biases in the seasonal cycle (see Sect. 3.2).
The observed synoptic variability is always larger than in the hindcast. By using monthly mean NEE the simulated variability is further dampened. This suggests that although the transport plays a first order role in the synoptic variability of atmospheric CO 2 , the day-to-day variability of NEE also plays an important role in enhancing it. This is confirmed in Table 4 where the correlations between the de-trended CO 2 from the model and observations at the two levels of the tall tower at Park Falls are shown for the two simulations with and without NEE day-to-day variability. The simulated CO 2 always correlates better with observations when the synoptic variability of NEE is included, except when the observations are sampling the free troposphere. That is the case for the 396 m level during nighttime, when large-scale advection dominates the variability and both simulations have very high correlations coefficients.
The passage of frontal low pressure systems is responsible for the long-range transport of CO 2 via their warm conveyor belts which lift CO 2 rich air from the surface to the mid and upper-troposphere. This large-scale advection is illustrated in Fig. 9 where positive CO 2 anomalies originating from the surface are shown in the region of frontal ascent at different vertical levels (850, 500 and 300 hPa). On 21 and 23-24 September, Park Falls experiences the advection of positive Atmos. Chem. Phys., 14, 11959-11983 CO 2 anomalies associated with the passage of two different low pressure systems. The cloudy warm conveyor belts in the mid-latitude low pressure systems are also associated with changes in temper-ature and solar radiation at the surface which in turn produce an increase in NEE (Fig. 7). This increase in NEE can be associated with a decrease in GPP following a decrease in radiation (e.g. 3 and 7 September), an increase in R eco follow- ing an increase in temperature (e.g. 21 September), or both a simultaneous decrease in GPP and increase in R eco due to a concurrent decrease in radiation and increase in temperature (e.g. 11 and 23-24 September). It is also interesting to note that on 29 September, the passage of a low pressure system lead to an increase in temperature at Park Falls, resulting in a simultaneous increase in GPP and R eco . In the model the increase in GPP is larger than the increase in R eco , leading to a decrease in NEE. This NEE decrease opposes the observed increase in atmospheric CO 2 .

CO 2 diurnal cycle
The diurnal cycle is assessed at two ICOS sites, one in Europe (Cabauw, the Netherlands) and one in Africa (Lamto, Ivory Coast). The amplitude of the diurnal cycle varies strongly at synoptic scales as shown by Figs. 10 and 12. This variability affects mainly the higher-values of CO 2 at nighttime, whereas the daytime CO 2 has a much lower monthly standard deviation (Fig. 11). As expected, the amplitude of the diurnal cycle decreases rapidly with height at the ICOS tall tower at Cabauw, Netherlands. The CO 2 hindcast is able to reproduce the changes in the amplitude of the diurnal cycle, both in time and in height. At the lower level (20 m), the model overestimates the variability of the nocturnal CO 2 values by largely overestimating the CO 2 peaks during three specific nights (24)(25)(26)(27). These are days when the 10 m wind speed drops to 1 m s −1 and the boundary layer height is very shallow. Under these conditions the CO 2 hindcast is highly uncertain because of both uncertainties in the mixing under stable conditions (Sandu et al., 2013) and the strong influence of the errors in the surface fluxes when the boundary layer collapses. In the hindcast, the daytime CO 2 trough is consistently underestimated at all vertical levels, which is consistent with the negative global bias described in Sects. 3.1 and 3.2.
At the tropical African site of Lamto, Ivory Coast, the diurnal cycle also shows the largest errors are at nighttime with an overestimation of CO 2 (Fig. 12), which is consistent with the positive bias in the CO 2 hindcast during the dry season. Nevertheless, it is clear that the nighttime overestimation does not occur every day (Fig. 12a). This suggests there Atmos. Chem. Phys., 14, 11959-11983  is a variable forcing responsible for the errors associated with the CO 2 hindcast. Correlations of the daily mean CO 2 with both boundary layer height and NEE fluxes from the model have been computed, in order to find which one is the main driver in the synoptic variability of the diurnal cycle amplitude. The daily mean boundary layer height from the model correlates well with the observed and modelled diurnal cycle amplitude of CO 2 at Cabauw with a correlation coefficient of −0.73 for the two of them. Both nighttime and daytime boundary layer heights play a role in the synoptic variability of diurnal cycle at Cabauw. At Lamto the most important factor explaining the synoptic variability of the diurnal cycle amplitude is the nighttime boundary layer height, with correlation values of −0.50 and −0.67 for the observed and modelled amplitude of the CO 2 diurnal cycle respectively. The correlation of the daily mean CO 2 and the NEE fluxes is below 0.3 at both sites. This implies that the NEE fluxes alone are not able to explain the synoptic variability of the diurnal cycle at those sites in September 2011. Although the boundary layer height at both Cabauw and Lamto appears to be the main factor explaining the variability of the diurnal cycle amplitude, this does not mean that the surface fluxes do not contribute. In fact, this evaluation shows that the surface fluxes and their errors have their effect enhanced under very stable conditions, when the boundary layer is very shallow.

Interhemispheric gradient of CO 2
The interhemispheric gradient is an important feature for CTM simulations, because it can be used to detect errors in both transport and CO 2 fluxes. As the TRANSCOM evaluation showed a good interhemispheric gradient for CH 4 in the IFS (P. Patra, personal communication, 2012), we expect most of the error to come from the CO 2 fluxes. The interhemispheric gradient of CO 2 has been evaluated using the HIPPO flight campaign data (Wofsy, 2011;Wofsy et al., 2012) in 2010 and 2011 (Fig. 13). In order to compare the simulated and observed CO 2 , the nearest model grid point, model level and model 3-hourly archived time to the observation is used. In March and April the comparison shows that the CO 2 -rich outflow from Asia in the region of the subtropical jet is overestimated in the simulations. Background biases fall between −1 and −4 ppm, except for the mid and high latitudes where the background biases range between  Figure 9. Transport of atmospheric CO 2 anomalies associated with the passage of low pressure systems over North America. The colours depict the CO 2 anomalies anomalies above the well-mixed background CO 2 at different vertical levels: grey near the surface, cyan at 850 hPa, blue at 500 hPa and dark grey at 300 hPa. The anomalies are defined as CO 2 dry molar fraction above the background value of 392 ppm for both near the surface and at the 850 hPa level; and above the background value of 388 ppm for the 500 and 300 hPa levels. The location of the TCCON sites are depicted by a red triangle. The black contours of mean sea level pressure show the location of the centre of the low pressure systems (L).
presented in Sect. 3.2. As a result of this negative bias in the lower mid-troposphere at NH mid-latitudes, the interhemispheric gradient is too strong in the summer and too weak in the spring. Similarly the negative vertical gradient between the lower and upper-troposphere in spring is too weak and the positive vertical gradient in the summer is too strong.

Vertical gradient of CO 2
One of the most important and more uncertain parts of the transport is the vertical mixing and the resulting vertical profiles over continental regions with strong surface fluxes (Kretschmer et al., 2012). There is a large variability between models in the simulation of vertical gradients and this strongly affects their consensus in the optimized NEE fluxes derived from different flux inversion systems (Stephens et al., 2007). In order to assess the performance of the hindcast in representing the vertical profiles, the model has been compared with observed vertical profiles at midday from NOAA/ESRL aircraft data in North America (Fig. 14a), following Stephens et al. (2007).
Results show an underestimation of the vertical gradient in both the lower and mid troposphere during winter ( Fig. 14b-d). The observed difference in the lower troposphere between altitudes of 1 and 4 km is +2.26 ppm compared to the modelled difference of +1.10 ppm. In the mid troposphere the discrepancy is smaller, +0.99 ppm between 4 and 6 km in the observations vs. +0.78 ppm in the model. The gradient is reversed and less steep during the summer. This is due to the change of sign in the NEE flux -from net release in winter to net uptake in summer -as well as the stronger vertical mixing associated with more convectively unstable atmospheric conditions. The model is able to simulate these changes, but still underestimates the observed gradient of −0.86 ppm in the lower troposphere compared to −0.47 ppm in the model.

Discussion
The hindcast performance is discussed in this section and possible ways of improving its deficiencies are described. The errors in the simulated CO 2 are dominated by errors in the fluxes. This is shown by the errors in the global budget, correlation coefficients and consistent biases computed us-Atmos. Chem. Phys., 14, 11959-11983 ing flight vertical profiles, total column observations as well as surface observations. The largest atmospheric CO 2 biases are in the NH, particularly in the Arctic region (north of 66 • N). However, this does not imply that the error in the fluxes is largest there. It is very likely that the larger negative biases in the arctic reflect the fact that the CO 2 biases from NH mid-latitudes (defined here between 30 • N and 66 • N) are transported northwards, consistently with the amplification of seasonal cycle in the arctic due to the coupling between mid-latitude fluxes and transport as described by Parazoo et al. (2011). The flux signal in the NH is coming predominantly from mid-latitudes, which include the boreal forests. Keppel-Aleks et al. (2011) demonstrated that small errors in NEE fluxes in the boreal region between 45 and 65 • N have a larger impact on the seasonal cycle amplitude of total column atmospheric CO 2 than changes at lower latitudes, due to the greater seasonality of NEE in the boreal region.
NH Spring is the season where the largest errors occur, both in budget (bias) and in the synoptic variability (correlations). Other models also found the spring months to have the lowest correlation coefficients with observed daily CO 2 (e.g. Geels et al., 2004;Pillai et al., 2011). This is not surprising as the onset of the CO 2 drawdown associated with the growing season causes a rapid shift in the dominant component of the NEE, i.e. from R eco to GPP. The simulated biogenic fluxes experience this shift in early spring (March) for two sites associated with cold temperate deciduous forest and corn crops; whereas in reality this shift occurs later on between April and May (see Fig. 8 of Falge et al., 2002). The intercomparison of several modelled NEE data sets with TCCON observations by Messerschmidt et al. (2013) showed that the best fit with the TCCON data was given by the SiB model which had the 20-75 • N aggregated NEE shift in April. The reasons for the one month error in the start of the growing season need to be further investigated. are the representation of the sensitivities of R eco and GPP to the variations of temperature and radiation in the CTESSEL model (Balzarolo et al., 2014), and the uncertainties associated with the estimation of the reference respiration as well as the simplistic radiative transfer scheme for vegetation (Boussetta et al., 2013a). Other seasons show much larger correlation coefficients, particularly in the NH winter and autumn where the variability is explained by the coupling between meteorology (i.e. transport) and flux variability associated with the passage of frontal low pressure systems. The correlation coefficients at the tall towers which are influenced by vegetation are higher than those presented so far by other models with similar or higher horizontal resolution (Geels et al., 2004;Pillai et al., 2011). This is very encouraging and it emphasizes the importance of the interaction between meteorological transport and forcing of the fluxes in the simulation of the CO 2 synop- tic variability. The NEE synoptic variability plays an important role, enhancing the CO 2 day-to-day variability locally. Within the boundary layer this effect is even more important, as the local fluxes play a more prominent role in modulating the atmospheric concentrations. In other words, the synoptic variability of atmospheric CO 2 could not be properly represented using climatological CO 2 fluxes, or offline biogenic CO 2 fluxes forced by climatologies of meteorological fields.
The sign of the vertical gradient is well represented in the hindcast, but the magnitude of the gradient is underestimated in the lower troposphere, particularly during winter. Although the fluxes can also be responsible for this underestimation, it is very likely that the vertical diffusion in the model is also contributing by having too much vertical mixing. This is a well-known problem in NWP models -including the IFS -which enhance the turbulent diffusion in stable conditions in order to compensate for errors caused by other poorly represented processes, such as orographic drag and the strength of the land-atmosphere coupling (Sandu et al., 2013).
The evaluation of the diurnal cycle also confirms that the boundary layer height and the 10 m wind speed are important controlling factors on the large day-to-day variability in the skill of the CO 2 hindcast. Under stable conditions when the boundary layer is shallower, there is an enhanced impact of the surface flux and their associated errors on the atmospheric CO 2 close to the surface. At the same time, the errors associated with turbulent mixing are also largest in stable conditions (Sandu et al., 2013).

Summary and further developments
This paper documents a new CO 2 forecast product from the MACC-II project, the pre-operational Copernicus atmospheric service. The CO 2 hindcast skill has been assessed at global to local scales and at temporal scales ranging from interannual variability to the diurnal cycle using a wide range of observations. Overall the hindcast can simulate very well the CO 2 synoptic variability modulated by the coupling between meteorological forcing of the fluxes and transport.
Comparing the synoptic variability with and without day-today variability in NEE indicates that in order to improve the synoptic skill of a CO 2 forecast, it is imperative to include and improve the day-to-day variability of the NEE fluxes, as well as its large-scale gradient. Improvements in the mod-elling of CO 2 fluxes and transport are expected as part of the ongoing efforts to upgrade the real-time CO 2 forecasting system of the Copernicus atmospheric service, in line with the updates of the operational IFS at ECMWF. For instance, the new developments in the convection and vertical diffusion parameterizations (Bechtold et al., 2014;Sandu et al., 2013) have been shown to have a positive impact on the diurnal cycle of convection and near-surface winds in the new IFS model cycle CY40R1. These improvements in the transport are also expected to lead to improvements in the CO 2 forecast. There are also developments in the assimilation of new satellite products in the IFS that could have a significant impact on the modelling of the CO 2 fluxes. For example, the assimilation of the near-real time albedo and LAI from the Copernicus Global Land Service , and   (Tans et al., 1996) in black and CO 2 hindcast in blue from 2003 to 2007 for January to March, January to December and July to September respectively. The average profiles plus/minus their standard deviation are shown as dashed lines. the SMOS/ASCAT soil moisture products (Muñoz-Sabater et al., 2012de Rosnay et al., 2012) could improve the phenology and the meteorological forcing on the modelled NEE fluxes, respectively. Further improvements of the vegetation radiative transfer scheme based on Carrer et al. (2013) are also planned for the near future.
Currently, the forecast is not constrained by CO 2 observations. Thus, there is an accumulating global bias (ranging from 2 to 4 ppm in magnitude). The bias is largest in the NH and it is associated predominantly with errors in the NEE fluxes in NH mid-latitudes, particularly during the growing season in spring. This model bias is larger than the bias of the currently available satellite CO 2 retrievals from GOSAT of only a few tenths of ppm (Notholt et al., 2013). Therefore, when such retrievals can be assimilated in near-real time in order to produce a CO 2 analysis to initialize the CO 2 forecast with, the bias of the forecast will also be reduced. Because the CO 2 forecast has good skill in simulating the synoptic variability of CO 2 in real time, it should provide a good background state for the assimilation of the available CO 2 observations and satellite retrievals from GOSAT, as well as other upcoming satellite missions, e.g. OCO-2.
The CO 2 observations provided in near-real time by the operational ICOS network are invaluable for the necessary CO 2 forecast error assessment. Continuous monitoring of the MACC-II CO 2 forecast based on the operational ICOS network is available online one day behind real time (www. copernicus-atmosphere.eu/d/services/gac/verif/ghg/icos).
The CO 2 forecast presented in this paper aims at providing information on the spatial and temporal variations of atmospheric CO 2 in real time. As such, it can be useful for a variety of purposes. For example, the atmospheric CO 2 fields can provide a link to collocate the CO 2 retrievals from satellite observations in time and spatial with ground-based observations for calibration, bias correction and evaluation purposes . Some satellite retrievals also rely on model-based CO 2 products to infer methane total columns and therefore avoid the expensive simulation of radiative scattering (Frankenberg et al., 2011). Other uses include the provision of boundary conditions for regional modelling and flux inversions (Matross et al., 2006;Rivier et al., Atmos. Chem. Phys., 14, 11959-11983, 2014 www.atmos-chem-phys.net/14/11959/2014/ 2010; Schuh et al., 2010;Broquet et al., 2011), helping the interpretation of observations (Schneising et al., 2012) and supporting the planning of field experiments (Carmichael et al., 2003). Finally, having real-time estimates for atmospheric CO 2 abundances has also other potential benefits, including a better representation of the model radiation and the radiance observation operator Engelen and Bauer, 2011), as well as evapotranspiration (Boussetta et al., 2013b) in NWP analyses and forecasts.
The real-time global CO 2 forecast is now part of the MACC-II suite of products freely available from the MACC-II data catalogue (http://www.copernicus-atmosphere.eu/ catalogue/). Before comparing the model with observations, the atmospheric CO 2 modelled fields need to be processed in order to match the modelled qualities with the observed quantities, including a collocation in space and time. The first step is to extract the vertical profile from the 3-hourly archived CO 2 forecast fields at the nearest land gridpoint to the location of the observation. For in situ observations, a linear interpolation to the observation height above the surface is performed in altitude. The last step is to collocate the observation and model in time. This is done by linearly interpolating the forecast data in time to match the observation time. The specific computations for the in situ (including the NOAA/ESRL flights) and total column observations are described below.

A1 In situ observations
In order to interpolate the model data to the sampling height for the in situ observations, the pressure of the model layer boundaries p l is converted at each grid point to altitude z h l by z h l = z h l+1 + R d g T l (1.0 + 0.61 q l ) ln p l+1 p l , where R d = 287.06, g = 9.8066 and l ranges from 1 to NLEV+1 (number of model levels+1) with z h NLEV+1 = 0. Then the elevation in the middle of the model layer z is computed for each level i ranging from 1 to NLEV by

A2 TCCON observations
The TCCON retrieved total columns are directly compared to the integrated averaging kernel-smoothed profile derived from the model CO 2 dry molar fraction profile (x m , directly extracted from the model without any correction required) following Rodgers and Connor (2003) and Wunch et al. (2010): where c s is the smoothed model forecast column average, c a is the a priori total column, a is a vector containing the TCCON absorber-weighted column averaging kernel, h T is a dry-pressure weighting function, and x a is the a priori CO 2 dry molar fraction profile. All the quantities of Eq. (A3) are interpolated onto the same model vertical grid. As the IFS has a hybrid-sigma pressure vertical grid, the model levels have corresponding pressure levels that vary in space and time.
The number of model levels (NLEV) used by the model forecast in this paper is 60. Note that the model does not provide any CO 2 dry molar fraction value at the surface. The model vertical levels are bounded by NLEV+1 half pressure levels (from 0 Pa to the surface pressure).

Equation (A3) can be re-written as
where c ak m and c ak a are the averaging kernel-weighted drypressure-weighted vertical columns from the model and a priori profiles respectively. The three terms are computed as the sum over each pressure level i: Note that h is an approximation of the dry-pressure weighted function following O'Dell et al. (2012) given by the following: with where q i and g i are the specific humidity and the gravitational acceleration at pressure level i (q i = q(p i ), g i = g(p i )), and M dry air is the molar mass of dry air. The bar denotes the average over a pressure layer computed as follows: