Greenhouse gas simulations with a coupled meteorological and transport model : the predictability of CO 2

A new model for greenhouse gas transport has been developed based on Environment and Climate Change Canada’s operational weather and environmental prediction models. When provided with realistic posterior fluxes for CO2, the CO2 simulations compare well to NOAA’s CarbonTracker fields and to near-surface continuous measurements, columns from the Total Carbon Column Observing Network (TCCON) and NOAA aircraft profiles. This coupled meteorological and tracer transport model is used to study the predictability of CO2. Predictability concerns the quantification of model forecast errors and thus of transport model errors. CO2 predictions are used to compute model– data mismatches when solving flux inversion problems and the quality of such predictions is a major concern. Here, the loss of meteorological predictability due to uncertain meteorological initial conditions is shown to impact CO2 predictability. The predictability of CO2 is shorter than that of the temperature field and increases near the surface and in the lower stratosphere. When broken down into spatial scales, CO2 predictability at the very largest scales is mainly due to surface fluxes but there is also some sensitivity to the land and ocean surface forcing of meteorological fields. The predictability due to the land and ocean surface is most evident in boreal summer when biospheric uptake produces large spatial gradients in the CO2 field. This is a newly identified source of uncertainty in CO2 predictions but it is expected to be much less significant than uncertainties in fluxes. However, it serves as an upper limit for the more important source of transport error and loss of predictability, which is due to uncertain meteorological analyses. By isolating this component of transport error, it is demonstrated that CO2 can only be defined on large spatial scales due to the presence of meteorological uncertainty. Thus, for a given model, there is a spatial scale below which fluxes cannot be inferred simply due to the fact that meteorological analyses are imperfect. These unresolved spatial scales correspond to small scales near the surface but increase with altitude. By isolating other components of transport error, the largest or limiting error can be identified. For example, a model error due to the lack of convective tracer transport was found to impact transport error on the very largest (wavenumbers less than 5) spatial scales. Thus for wavenumbers greater than 5, transport model error due to meteorological analysis uncertainty is more important for our model than the lack of convective tracer transport.


Introduction
Atmospheric observations of CO 2 are important for understanding the global carbon cycle and its response to perturbations from anthropogenic emissions into the atmosphere.The global carbon budget is routinely updated using atmospheric and oceanic measurements in conjunction with a careful accounting of anthropogenic emissions (Le Quéré et al., 2015).Since the greatest uncertainty in the quantified exchanges of carbon between the atmosphere, ocean and land reservoirs is associated with the terrestrial biosphere flux (Friedlingstein et al., 2014), this component is determined as the residual between anthropogenic emissions and atmospheric and ocean changes (which are relatively well constrained by measurements).The growth in global atmospheric CO 2 is attributed to anthropogenic emissions but superimposed on this trend are seasonal and interannual variations which are largely attributed to the terrestrial biosphere (Ciais et al., 2013).Specifically, observations of the 13 C isotope (which is fractionated during photosynthesis by land plants so its relative concentration is indicative of the terrestrial biosphere) reveal that the interannual and seasonal variations of the global atmospheric CO 2 budget are due to the terrestrial biosphere (Keeling et al., 2005).At the same time, biospheric uptake is influenced by climate variations (Nemani et al., 2003;Friedlingstein et al., 2006;Zhao et al., 2011).Specifically, the El Niño-Southern Oscillation (ENSO) signal and volcanic eruptions can explain 75 % of the variations in the terrestrial biospheric uptake (Raupach et al., 2008).Thus, atmospheric observations are important for determining the global terrestrial biospheric flux, while temporal biospheric flux variations explain atmospheric variability of CO 2 concentrations on seasonal and interannual timescales.
Beyond the global budget, an important challenge is to understand how temporal variations in atmospheric sources and sinks reflect the interplay between natural processes and anthropogenic perturbations and also the feedback between the carbon cycle and climate variations.By combining atmospheric observations with atmospheric models, a spatial distribution of fluxes of CO 2 between the Earth's surface and the atmosphere can be obtained through inverse modelling (e.g.Rödenbeck et al., 2003;Patra et al., 2005;Baker et al., 2006a;Peylin et al., 2013;Chevallier et al., 2014).Such "flux inversions" performed on the global domain and incorporating only around 100 or so CO 2 observation stations near the surface are able to constrain the global atmospheric carbon budget, capture interannual and seasonal variations and attribute these to the terrestrial biosphere (Rödenbeck et al., 2003;Peylin et al., 2013).In fact, only a few observation sites representative of background values (far from CO 2 sources and sinks) are needed to constrain the global CO 2 budget because variations in the vicinity of source/sink regions are smoothed out by the time these locations are reached (Bruhwiler et al., 2011;Keeling et al., 2005).However, with more surface observations, the retrieved flux uncertainties can be reduced and the ability to retrieve smaller spatial scales is improved (Bruhwiler et al., 2011).Moreover, with the desire to understand the interplay between the natural and anthropogenically perturbed processes, observations near source and sink regions become important.Such observations will be influenced by atmospheric variations on diurnal, synoptic, seasonal and interannual timescales.On the diurnal timescale, the variation of CO 2 due to uptake by plants through photosynthesis in sunlit hours is strongly modulated by turbulent transport through the planetary boundary layer (PBL), which also evolves throughout the day.This is the so-called rectifier effect which helps to explain the annual mean north-south gradient of CO 2 (Denning et al., 1995).Specifically, the uptake of CO 2 by plants during the spring growing season occurs when the PBL is generally unstable and deeper, while in winter CO 2 increases due to biospheric respiration can build up when the boundary layer is stable and shallow.In addition, poleward heat transport by baroclinic disturbances is stronger in winter in the northern extratropics, preferentially transporting high CO 2 values north relative to the summer (Chan et al., 2008).Synoptic-scale systems also influence CO 2 evolution, particularly in northern midlatitudes where advection can explain up to 70 % of day-to-day variability (Parazoo et al., 2008).On the interannual timescale, variations in biospheric uptake can be partially attributed to climate variations (Patra et al., 2005) and flux inversion systems are able to attribute interannual variability of the global CO 2 budget to the tropical biosphere (Baker et al., 2006a).Because atmospheric observations of CO 2 contain the signals of both surface fluxes and atmospheric variations, and they are needed for data assimilation in state estimations or flux inversions, it is important to be able to accurately characterize and model how the atmosphere modulates CO 2 evolution.
In a flux inversion system, observations of CO 2 concentration are used to solve for surface fluxes.To relate the surface flux to an atmospheric concentration, an offline atmospheric transport model is typically used.In doing so, the mismatch between observed and modelled concentrations can be inverted to estimate a flux increment, but the atmospheric model's transport is assumed to be perfect (Baker et al., 2006b).The fact that it is not means that one source of uncertainty in the inverse problem is not accounted for.For this reason, the need to characterize "transport errors" has been well recognized and has led to the formation in 1993 of an international group called TransCom (http://transcom.project.asu.edu/),focused on understanding and quantifying the contribution of transport errors to flux estimates from inverse models.Because different models are used by different flux inversions, retrieved fluxes based on multiple inversions may be more reliable than those from a given system since the unknown transport error (if it is random) may be averaged out.Moreover, the sensitivity of flux inversion results to transport error can be identified since different models and hence different transport errors are used.Transport error generally refers to the deviation of CO 2 prediction from an (un-knowable) true value due to the use of an imperfect transport model and thus includes errors due to model formulation (e.g.convective or boundary layer parameterizations or advection schemes), the use of imperfect meteorological analyses and uncertain CO 2 initial conditions.Transport errors have been found to be an important source of errors in flux inversions (e.g.Chevallier et al., 2014Chevallier et al., , 2010;;Houweling et al., 2010;Law et al., 1996).The sources of transport error arising from model formulation include errors in the representation of mixing in the planetary boundary layer (Denning et al., 1995), vertical mixing in the free atmosphere (Stephens et al., 2007;Yang et al., 2007), synoptic-scale and frontal motions (Parazoo et al., 2008) and convective transport (Ott et al., 2011;Parazoo et al., 2008).
The uncertainty in meteorological analyses is also an important source of transport error (Liu et al., 2011).In an effort to address this type of uncertainty which cannot be accounted for in flux inversion systems, NOAA's CarbonTracker attempted to perform an ensemble of inversions using different atmospheric analyses as well as different prior flux sets. 1ith a coupled meteorological and tracer forecast model (as used by Liu et al., 2011), the impact of meteorological uncertainties on tracer transport is more easily identified.Recently, coupled meteorological/tracer forecast models have been developed at operational centres such as the European Centre for Medium Range Weather Forecasting (ECMWF) (Agustí-Panareda et al., 2014) and NASA Goddard's Global Modelling and Assimilation Office (GMAO) (Ott et al., 2015).In these operational systems, short-term predictions of greenhouse gases are produced and satellite observations are assimilated at ECMWF (Massart et al., 2016).Such products may provide useful background or a priori information for satellite retrievals as well as providing initial and boundary conditions for regional flux inversions.Coupled meteorological and CO 2 data assimilation systems also provide useful information on the reliability of correlations between meteorological and CO 2 transport errors (Kang et al., 2011) and on the temporal propagation of the observed signal in the context of transport errors (Kang et al., 2012).
The goal of this work is to better understand the component of transport error that is due to uncertain meteorological states.While Liu et al. (2011) have shown that meteorological forecast errors in the presence of CO 2 gradients produce CO 2 transport errors, here we consider the spatial scales identifiable in the context of imperfect meteorological analyses.Predictability concerns the study of forecast uncertainty and thus corresponds to the study of transport model errors.While transport errors comprise flux errors, model formulation errors, initial state errors and meteorological state errors, it is the latter that is the focus of this work.Using a coupled meteorological and tracer transport model, we first study the loss of CO 2 predictability due to uncertain meteorological initial conditions on weather (2 weeks) and seasonal timescales in order to obtain an upper limit to forecast errors arising from meteorological analysis errors.Then the spatial scales of errors in CO 2 arising from the use of imperfect atmospheric analyses are determined.It is shown that there is a spatial scale below which CO 2 cannot be retrieved simply due to the presence of meteorological analysis uncertainty.The advantage of isolating and comparing different components of transport error for a given model is then demonstrated.The spatial scales of model errors (such as the lack of convective transport of tracers) are compared to the scales retrievable in the context of imperfect atmospheric analyses.
The article is organized as follows.In Sect.2, we discuss the concept of predictability for the tracer transport problem and identify the components of transport error, while in Sect. 3 we describe a new model for greenhouse gas transport based on an operational weather and environmental prediction model.In Sect.4, this model is assessed in terms of its meteorology (24 h forecasts are compared to reanalyses used in flux inversions) and its CO 2 transport.The component of CO 2 transport model error due to meteorological initial conditions is studied first in Sect.5, considering both weather and seasonal timescales.This is followed by a study of the spatial scales of transport model error arising from uncertain meteorological analyses.Results are discussed and summarized in Sect.6.

Predictability and transport error
Atmospheric predictability concerns the quantification of forecast uncertainty with a weather or climate prediction model.The limited ability to forecast weather is related to the underlying chaotic dynamics of the atmosphere and is evident in the sensitivity of forecasts to initial conditions.However, even though the predictability of weather does not extend beyond 2 weeks, components of the climate can be predicted on seasonal and interannual timescales.The sources of predictability on climate timescales are related to atmospheric boundary conditions such as sea surface temperature, soil moisture, snow cover, vegetation and sea ice (Shukla and Kinter III, 2006).Because CO 2 is transported by the atmosphere, the predictability of the atmosphere is directly relevant to the predictability of CO 2 .However, atmospheric CO 2 concentrations are also strongly influenced by surface sources and sinks.In fact, when concentrations are globally averaged to compute a global CO 2 budget, it is only these fluxes that are relevant since atmospheric transport does not change global mass.However, when considering the spatial distribution of CO 2 , both surface fluxes and atmosphere transport contribute to the quality of predictions of CO 2 .While the forecast of CO 2 is not of such great public inter-S.M. Polavarapu et al.: Greenhouse gas simulations: the predictability of CO 2 est as a weather forecast, the ability to predict CO 2 is needed for solving flux inversion problems.Specifically, a CO 2 prediction is compared to atmospheric observations to retrieve fluxes.Thus knowledge of the accuracy of the CO 2 prediction is needed in the formulation of the model-data mismatch error statistics needed for solving the flux inversion problem (see Baker et al., 2006b).In addition, once the fluxes have been inverted, the transport model is integrated and compared to measurements to assess the quality of the flux inversion.Thus the quality of the CO 2 prediction is assessed.
Since CO 2 prediction is frequently performed with a transport model, particularly in the context of flux inversions, the prediction or forecast error is called "transport error".There is considerable concern over the nature and quantification of various sources of transport error (see http://transcom.project.asu.edu/).As noted in the Introduction, atmospheric variability on diurnal, seasonal and interannual timescales impact CO 2 concentrations.Thus atmospheric predictability on all of these timescales impact the predictability of CO 2 .Indeed, transport errors on diurnal timescales due to mixing in the boundary layer (Denning et al., 1995), synoptic timescales (Chan et al., 2004;Parazoo et al., 2008), seasonal timescales (Denning et al., 1995;Chan et al., 2008) and interannual timescales (Baker et al., 2006a) have all been identified as contributing to flux estimation errors.
The error of a CO 2 transport model can be defined as the departure of a CO 2 prediction from the true (but not known) value.Mathematically, this is expressed as , where the updated concentration at time step k+1 is given by a transport model: where T k,k+1 is the transport model which evolves the meteorological state, x k , the concentrations, c k , and the flux, s k at time step k to time step k + 1.This transport error has contributions from various sources which can be identified by considering the hypothetical case of a forecast obtained with perfect knowledge of fluxes, meteorology and initial CO 2 state.Such a forecast would still be imperfect because the transport model is imperfect: The last term is the model error which arises due to errors in the model formulation.Various authors have investigated the role of model errors due to boundary layer mixing (Denning et al., 1995), vertical mixing in the free atmosphere (Stephens et al., 2007;Yang et al., 2007), synoptic and frontal motions (Parazoo et al., 2008), convective transport (Ott et al., 2011;Parazoo et al., 2008) and mass conservation schemes (Agusti-Panareda et al., 2016).Note that all but the last of these model errors also impact meteorological forecast uncertainty in the case of a coupled meteorological and tracer transport model.Of course, model error is only one component of transport error.The additional error sources are due to errors in the meteorological state ε (1) and ( 2) into the expression for transport error and using these additional error definitions reveals that transport error is the convolution of all of these error sources.In other words, where Thus, the various error subcomponents interact with each other to produce transport error.When an offline transport model is used in a flux inversion, the meteorological state is assumed to be known so the first and last terms in Eq. ( 3) are absent.Since the meteorological state is obtained from a meteorological analysis (usually from a foreign assimilation or reanalysis system) the meteorological state error evolves from an analysis error.It is difficult to account for meteorological state errors on CO 2 transport errors in the flux inversion problem aside from producing an ensemble of inverse results with different meteorological forcings (e.g. from different reanalyses).However, with a coupled meteorological and transport model, Liu et al. (2011) have demonstrated that meteorological uncertainty is an important component of transport error (1.2-3.5 ppm at the surface and 0.8-1.8ppm in a column mean).Since the goal of flux inversion is to attribute the mismatch of CO 2 predictions and observations to a flux adjustment, a good understanding of the nature and size of all of the components of transport error is needed because only then can a successful attribution of the mismatch to flux be made.While progress has been made on identifying the role of model and flux errors on transport errors, there is less known about the impact of meteorological state errors on CO 2 predictability.With a coupled meteorological and transport model, we are able to consider the impact of meteorological state errors on CO 2 predictability and thus on CO 2 transport error.

Model description
In order to simulate greenhouse gas evolution, we can take advantage of the comprehensive forecast models already developed for operational environmental prediction at Environment and Climate Change Canada (ECCC).However, it was necessary to adapt these models for the purpose of tracer transport on multiannual timescales.Specifically, it was necessary to implement a mass conservation scheme, redesign the tracer variables, modify the vertical mixing in the boundary layer and add convective tracer transport.The basic modelling tools are described in the next subsection while the adaptations needed for this work are presented in Sects.3.2 (mass conservation and tracer variable definitions), 3.3 (horizontal diffusion), 3.4 (convective transport) and 3.5 (boundary layer mixing).The coupled meteorology and tracer transport forecasting system is described in Sect.3.6.

The Canadian operational environmental prediction models
For many decades, the Canadian Meteorological Centre (CMC) has been producing operational weather forecasts for public dissemination.Since 24 February 1997, these forecasts have utilized the Global Environmental Multiscale (GEM) model (Côté et al., 1998a, b;Girard et al., 2014).GEM is a grid-point model which solves the hydrostatic (global domain) or nonhydrostatic (Yeh et al., 2002) (regional domain) primitive equations using a hybrid terrainfollowing vertical coordinate (Girard et al., 2014).As of February 2013, the grid spacing of the global model is roughly 25 km, originally using a regular lat-long grid and, since December 2015, a yin-yang grid (Qaddouri and Lee, 2011).There are 80 vertical levels spanning the surface to 0.1 hPa.The usual physical processes of radiation (Li and Barker, 2005), boundary layer mixing (Bélair et al., 1999), shallow (Bélair et al., 2005) and deep convection (Kain and Fritsch, 1990;Kain, 2004), orographic gravity wave drag (McFarlane, 1987) and nonorographic gravity wave drag (Hines, 1997a, b) are included in all model configurations.The land surface model and assimilation scheme are described in Bélair et al. (2003a, b).More details of the physics package are found in Mailhot et al. (1998).
Operational air quality forecasts have been produced by CMC since 2001 in order to provide real-time forecasts of the air quality health index on a limited area domain covering most of North America.As of 18 November 2009, these forecasts have utilized GEM-MACH (Modelling Air quality and CHemistry) (Moran et al., 2010;Robichaud and Ménard, 2014;Makar et al., 2015).GEM-MACH is a version of GEM in which complete tropospheric chemistry (involving over 100 chemical reactions) is modelled online, where "online" refers to the fact that the chemistry module is fully integrated into the meteorological model time step.The operational products involve an analysis of ground-level ozone, fine and coarse particulate matter (PM 2.5 and PM 10 ), NO 2 and SO 2 over a limited area domain covering North America (Robichaud and Ménard, 2014).The grid spacing currently used is 10 km horizontally, with 80 vertical levels from the surface to 0.1 hPa.The operational forecasts are driven by time evolving meteorological boundary conditions from the operational regional deterministic prediction system (Fillion et al., 2010;Caron et al., 2015), while the chemical boundary conditions are defined using predetermined seasonally-averaged states.A global version of GEM-MACH is also in development for the purpose of providing boundary conditions for the regional model and a parameterized stratospheric chemistry model (McLinden et al., 2000) is used for UV index forecasting.
Our primary interest is in global greenhouse gas distributions, and thus we chose to use the global GEM-MACH configuration.In the future, extensions of our system for regional greenhouse gas simulations could be based on the operational regional GEM-MACH configuration.By definition, an operational forecasting system is constantly changing.Since our period of interest commences in 2009 with the launch of GOSAT on 23 January and since the upper boundary of GEM was raised from 10 to 0.1 hPa on 22 June 2009 (Charron et al., 2012), this stratospheric model configuration was chosen.However, for computational expediency (especially during the model development phase) the grid spacing was coarsened to 0.9 • (400 × 200 grid points), which is roughly twice the grid spacing used by CMC's global deterministic prediction system (800 × 600 grid points) in 2009.The model time step is 15 min.
Our primary focus in this work is global carbon dioxide simulation, but extensions for other greenhouse gases (methane and CO) are also being developed.For that we employ a simple parameterized climate chemistry involving a single OH reaction for the methane and CO simulations.This chemical module is activated from the GEM-MACH chemical interface, which is also used to handle emission inputs.This simplified model version, which only includes the treatment of greenhouse gas (chemistry and transport), will be called GEM-MACH-GHG hereafter.

Mass conservation and tracer variable definitions
The dry air mass of the atmosphere is known to be constant, since changes in trace gases are very small (Trenberth, 1981;Trenberth and Smith, 2005).However, as with many weather and climate prediction models, GEM does not conserve dry air mass.Specifically, GEM loses 0.1 hPa in global mean surface pressure during a 10-day forecast, which is precisely the same rate as that seen in the ECMWF model (Diamantakis and Flemming, 2014).This loss is only 0.01 % of the global mass and thus there is a negligible impact on medium-range weather forecasts.However, for longer simulations relevant for climate timescales, this error can accumulate.Thus GEM (as does the ECMWF and other models) has a parameter to allow the global mean surface pressure to be conserved by adding a spatially uniform adjustment to each grid cell.The constant is determined by the constraint that the air mass at the end of the dynamics step equal that at the start of the step.For our simulations of tracer transport, it is necessary to use this switch to enforce the conservation of dry air mass.
The advection scheme in GEM uses a semi-Lagrangian approach as it affords longer time steps and the computational efficiency desirable in an operational context.How-ever, semi-Lagrangian schemes are well known to be nonconservative (Williamson, 1990;Staniforth and Côté, 1991).GEM is a grid point model, and the semi-Lagrangian approach involves first determining the upstream positions of the grid cells (using a 3-D trajectory in our case) and then interpolating the value of the advected field at these locations.During this first step, gradients in the wind field will limit the accuracy of advection and during the second step, errors associated with the interpolation of the advected field can deteriorate species mass conservation.By combining a finite volume approach with the semi-Lagrangian scheme, it is possible to devise inherently mass-conserving schemes such as SLICE (Semi-Lagrangian Inherently Conserving and Efficient) (Zerroukat and Allen, 2012).Indeed, this scheme is being investigated for implementation with GEM in the future.In the meantime, an interim solution is needed.The approach taken here, as in Diamantakis and Flemming (2014), is to adopt a global mass fixer.To be clear, the global mass is obtained by summing the individual masses of carbon in each grid box.The mass of carbon in a grid box is determined from the mixing ratio (mass of carbon divided by the mass of dry air) times the mass of dry air in the same box.A global mass fixer scheme computes the global mass of tracer at the beginning and end of the advection step and then distributes the change in global mass among various grid cells.There is no unique way of distributing the mass spatially, so different choices result in different mass fixer schemes.While a few approaches were tried and tested with a few different chemical species (ozone, CH 4 and CO 2 ), the one which produced the most physically desirable results was that by Bermejo and Conde (2002).In this scheme, the global mass is distributed according to the smoothness of the field.That is, mass is preferentially adjusted where gradients are larger (and interpolation error is known to be larger).The exact scheme used in GEM is precisely that selected by ECMWF and it is described in Diamantakis and Flemming (2014, Sect. 3.1).For greenhouse gas transport, the ECMWF forecast model also uses the Bermejo-Conde scheme (Agusti-Panareda et al., 2016).
For tracers, another issue with the interpolation step of the semi-Lagrangian advection scheme is the potential creation of spurious subgrid-scale structure due to Gibbs effects.To avoid this problem, GEM uses a quasi-monotonic interpolation scheme (Bermejo and Staniforth, 1992) which combines high-order and low-order interpolation schemes to prevent overshoots and undershoots, thereby preventing the formation of spurious extrema.However, just as Flemming and Huijnen (2011) noted, enforcing monotonicity was found to worsen the problem of mass non-conservation because it tends to diffuse gradients.Therefore, monotonicity was replaced by an iterative locally mass conserving (ILMC) scheme (Sørensen et al., 2013).The idea behind the ILMC scheme is to locally preserve the shape of the field and its gradients by distributing the excess (or deficit of) mass due to spurious extrema in the cubic interpolated field to ever in-creasing shells around the upstream departure point but in such a way that spurious extrema are avoided.By design, the use of ILMC does not impact our lack of global mass conservation, so it is used for preserving positive definiteness of tracer fields in place of a quasi-monotonic interpolation scheme.The ILMC is performed first, and then the Bermejo and Conde (2002) scheme is applied to the ILMC-corrected field.Details on the implementation of monotonic and massconservation schemes in GEM and their impact on species transport are found in de Grandpré et al. (2016).
The primitive equations solved by GEM are naturally written in terms of a moist density or pressure variable.However, when tracers are defined as mixing ratios with respect to moist air, tracers become coupled to the water vapour evolution.Thus one of two options must be taken.Either the tracer variable must be redefined every time the water vapour field is modified in the model or the tracers must be defined as mixing ratios with respect to dry air (e.g.Neale et al., 2010).The latter option has the advantage that measurements of CO 2 are frequently made in terms of dry mole fraction and are thus more easily related to a mixing ratio with respect to dry air.Since the nominal tracer equation within GEM assumed a moist mixing ratio, for this work GEM was modified to permit tracers to be defined as mixing ratios with respect to dry air (DMR).This involved (1) modifying the global mean surface pressure adjustment to ensure global dry air mass conservation, (2) converting the mass fixer scheme (i.e. the Bermejo-Conde scheme) and the ILMC to deal with DMR, (3) modifying the vertical diffusion equation to handle DMR and (4) ensuring emissions are correctly inserted into GEM's bottom model layer.Once all of these changes had been made, the tracer mass change in a given time step was found to still display a coupling to the water vapour change in that same time step.The reason was found to be due to the fact that the continuity equation in GEM does not account for the change in surface pressure due to the change in water vapour which occurs as a result of mass flux at the Earth's surface (i.e. because of precipitation or evaporation).By adding a new adjustment to surface pressure field after the physics step to account for the change in global water vapour mass, and correspondingly redefining the tracer variable using a vertical regridding approach similar to Jöckel et al. (2001) except for the use of dry pressure coordinates, global mass conservation of CO 2 was obtained during a model forecast.This deficiency in GEM's continuity equation is not unusual as it was also present in NASA GMAO's GEOS5 model (Takacs et al., 2015).The adjustment to the model dynamics done here to account for the modified continuity equation is rather similar to that described in Sect.2a of Takacs et al. (2015) although our work was done independently.The main difference is that only water vapour mass is considered here whereas cloud liquid water and cloud ice mass are also considered in Takacs et al. (2015) (although the authors note that additional masses due to liquid and ice phases are negligible).As these authors demonstrate, without the correction to the model's continu-ity equation, dry air mass conservation is not assured.They also note that in the context of an assimilation cycle, further constraints on dry air mass are needed during the analysis and initialization steps.Furthermore, they show that the impact of the errors is evident in the mismatch of global water vapour tendency based on the water vapour field itself vs. that based on precipitation minus evaporation changes for reanalysis products such as MERRA (Rienecker et al., 2011) and ERA-Interim (Dee et al., 2011).

Horizontal diffusion
It is standard to apply an explicit high-order (typically ∇ 6 ) diffusion operator to the meteorological fields in applications involving the global GEM model.It is arguably consistent therefore to apply the same level of diffusion to the tracer fields.However, because the constituent variables are defined as mixing ratios, the operator is not mass conservative.Because of the effort involved in obtaining a careful accounting of tracer mass (described in Sect.3.2), it was decided to not apply any horizontal diffusion to the tracers.

Convective transport
The parameterization of deep convection most frequently used by GEM and GEM-MACH is due to Kain and Fritsch (1990).The Kain and Fritsch (KF) scheme is based on a single column bulk mass flux approach.Entrainment of ambient air into the cloud environment associated with updrafts and downdrafts is proportional to the corresponding mass flux and inversely proportional to the cloud radius.So narrower convective towers experience more entrainment of lower-buoyancy ambient air and consequently have less intensity and lower cloud tops.The KF scheme is used in several forecast models (e.g.Japanese Meteorological Agency model, Saito, 2012; Bologna Limited-Area Model (BO-LAM), Lagouvardos et al., 2003;High-Resolution Limited-Area Model (HIRLAM), Eerola, 2013;Weather Research and Forecasting Nonhydrostatic Mesoscale Model (WRF-NMM), Gallus Jr. and Bresch, 2006).With the KF scheme, GEM has a good representation of convectively coupled waves and is able to capture the Madden-Julian oscillation (Lin et al., 2008).
The version of the KF scheme used for global deterministic prediction, which was our starting point, had to be modified for the purpose of greenhouse gas transport.The original parameters resulted in too frequent penetration of convection into the stratosphere in the tropics.This is because the updraft core radius had been set to 1500 km globally.This is a reasonable value for extratropical convection over land, and since the target region of ECCC's forecasting systems is Canada, it is a valid choice.However, observations (Lucas et al., 1994, and references therein) show that the updraft core radius varies with latitude and from land to ocean.Thus over the oceans, a value of 900 km is used in the trop-ics with 1000 km in the extratropics.Over tropical land grid points, a value of 1200 km is used.With these settings, tropical convection overshoots fall in the correct range (below 80 hPa).However, the model climatology is altered.The parameter changes impact land-sea contrasts which influence stationary Rossby wave generation and the spin-up of the Brewer-Dobson circulation in the stratosphere.They also affect flow over topography, which impacts orographic wave drag and its associated transport circulations.An animation (Fig. S1 in the Supplement) shows a comparison of the impact of changing KF parameters vs. that of adding convective tracer transport on CO 2 evolution.As expected, the change in parameters impacts the stratospheric distribution after a few months of simulation, whereas the introduction of tracer transport through deep convection impacts the tropical tropospheric CO 2 distribution at all times as well as the Northern Hemisphere in summer.The zonal mean values are small but so too are zonal standard deviations (both are less than 0.3 ppm).Another animation which compares the impact of changing KF parameters to that of adding convective tracer transport on column mean or XCO 2 (Fig. S2 in the Supplement) reveals that the magnitude of the impact is smaller for the change in parameters (maxima of 0.3 ppm) than for the introduction of tracer transport through deep convection (maxima of 0.8 ppm).

Vertical mixing
Turbulence in the PBL is important for the transport of heat, momentum, moisture and constituent fluxes from the surface to the atmosphere.The PBL scheme in GEM is described in Bélair et al. (1999), Benoit et al. (1989) and Mailhot and Benoit (1982).Since a summary of GEM's PBL scheme was recently presented in McTaggart-Cowan and Zadra (2015) and Aliabadi et al. (2016), we note here only the main differences between GEM's PBL parameterization and those of other models used for greenhouse transport or flux inversion.
Shallow convection also exerts a significant impact on the tracer distribution and vertical mixing in the lower troposphere (e.g.Dacre et al., 2007).The shallow convection scheme in GEM is a modified Kuo parameterization (Kuo, 1974) developed at ECCC (Bélair et al., 2005).As with the original Kuo scheme it does not address convective tracer transport.However, unlike the KF scheme it is not based on a mass flux formulation that can be easily adapted to handle tracers.For this reason the version of GEM used for this study does not include tracer transport by shallow convection.
The vertical transport of greenhouse gas fluxes follows that of other constituents and is described by where C is the constituent mixing ratio on resolved scales, ρ is air density, w is subgrid-scale vertical velocity and the product with an overbar is the vertical flux of constituent due to subgrid-scale turbulence.The first equality denotes the impact of subgrid-scale motions on the resolved constituent distribution.The subgrid-scale flux is parameterized through the second equality and is a function of the vertical gradient of the constituent.γ is a counter gradient flux relevant for an unstable PBL.K is the thermal eddy diffusivity and it depends on properties of the flow.Most models define K diagnostically in terms of an eddy length scale and local gradients of wind and virtual potential temperature.Such schemes are said to invoke a first-order closure.However, GEM's scheme specifies K through a prognostic equation for turbulent kinetic energy and is thus deemed a 1.5-order closure scheme (Holtslag, 2015).When K is determined diagnostically through local flow properties, the case of an unstable or strongly convective PBL is not well represented.Thus, an extension for non-local mixing due to large-scale eddies is often invoked in which the counter-gradient term represents large eddy flux and K permits eddy length scales comparable to the PBL height in the case of unstable boundary layers.CarbonTracker (Peters et al., 2004) and GEOS-Chem (Lin and McElroy, 2010) both use such a nonlocal scheme based on Holtslag and Boville (1993).ECMWF also uses a diagnostic nonlocal scheme (Köhler et al., 2011).In GEM's 1.5-order closure scheme, K can also represent nonlocal effects due to large eddies in an unstable PBL through a careful definition of eddy length scale (Bélair et al., 1999).One notable difference between GEM's prognostic nonlocal scheme and diagnostic nonlocal schemes is that the latter require PBL height as an input for determining K, whereas in GEM's scheme the PBL height is diagnosed from the turbulent kinetic energy profile (so K is not a direct function of PBL height).
A minimum value for eddy diffusivity of 0.1 m 2 s −1 is imposed in the PBL in GEM-MACH for air quality applications.This value was empirically chosen to balance the impacts on the various tropospheric reactive species (P.Makar, personal communication, 2012).In this work, the minimum value was raised to 10 m 2 s −1 because a value of 1 or 0.1 was found to occasionally lead to too little mixing.This was most evident in spuriously low CO 2 concentrations in the daytime during summer when biospheric fluxes are large and negative.(Note that the eddy diffusivity is kept fully varying in space -it is only the minimum value that is slightly altered.This is in contrast to earlier versions of GEOS-Chem in which instantaneous full mixing of emissions and mixing ratios occurred in the PBL.See Lin and McElroy, 2010, for example.)A negative consequence of raising the minimum value to 10 is a reduced amplitude of the diurnal cycle of CO 2 (as will be seen in Sect.4).However, many models have difficulty in capturing the amplitude of the diurnal cycle (Law et al., 2008;Patra et al., 2008) so this consequence was considered tolerable at present.Moreover, the comparison of model results to measurements on sub-diurnal timescales is very difficult, leading Law et al. (2008) to recommend that "com-parisons with observations should only be made for daily or longer time averages, and possibly for only part of the diurnal cycle".Thus, many inversion systems assimilate only afternoon mean observations and do not attempt to capture sub-diurnal timescales (Peters et al., 2010).

The coupled meteorology and tracer transport forecast cycle
To transport greenhouse gases, a simulation cycle is used.
Meteorological analyses archived at the CMC are inserted periodically to constrain the transport to reality.While operational assimilation cycles use an update frequency of 6 h, model forecasts during such a short period will be contaminated by spurious gravity wave generation, as evident in surface pressure time series (as well as other fields) (see Daley, 1991, chap. 6).However, after 24 h, spurious gravity waves have generally dispersed so an update frequency of 24 h was chosen, just as in Agustí-Panareda et al. ( 2014).(Note that a 6 h update cycle had also been tested and the resulting CO 2 evolution was very close to that obtained with a 24 h update cycle, but because the stratospheric circulation was poorer the 24 h update cycle is preferred.)A schematic diagram of the transport cycle is depicted in Fig. 1.The upper half of the figure depicts the operational system which collects meteorological observations over a 6 h window and uses these in a 4-D variational assimilation or 4D-Var (Gauthier et al., 2007) to generate an estimate of the meteorological state for a deterministic prediction.Although the operational deterministic prediction system now uses a hybrid approach to background error covariance estimation (Buehner et al., 2015), the analyses used here (from 2009 to 2010) were generated using the previously operational 4D-Var system (Charron et al., 2012).On the cycle's start date at 00:00 UTC, the meteorological analysis is combined with an initial state for greenhouse gases (here CO 2 only) and a 24 h coupled meteorology and tracer forecast is produced.The 24 h tracer forecast is subsequently combined with the meteorological analysis for the next day at 00:00 UTC (blue boxes) to produce a new coupled initial state for the second day's forecast (red filled circles).Note that no additional "initialization" or filtering (see Daley, 1991, chap. 6, 9, 10) scheme such as a digital filter or incremental analysis updates (Bloom et al., 1996) is used, as the impact on the CO 2 field was found negligible.Also note that in our simulation cycles, no assimilation of greenhouse gases is performed.An assimilation system is currently in development as discussed in Sect.6.In Fig. 1, every time a new analysis is inserted, a new surface pressure field informed by atmospheric observations is introduced.This surface pressure analysis will differ from the 24 h forecast of surface pressure because the model is not perfect, so the model forecast cycle will experience an abrupt shift in global mean surface pressure and, hence, global air mass.The change in global air mass will then impact the global tracer mass.To maintain global tracer mass conserva- Meteorological analyses were precomputed by CMC's operational global deterministic prediction system using all observations collected in a 6 h window centred on the analysis time.These analyses are represented by blue boxes.CO 2 tracers are added to these analyses to form an initial condition (solid red circles) for launching a 24 h forecast (red arrows) using the coupled meteorological/tracer model starting at 00:00 UTC of each day.tion across this discontinuity, it is necessary to redefine the tracer mixing ratio for the change in air mass.The scheme adopted follows that used by GEM for the global surface pressure adjustment.Specifically, a spatially uniform increment in tracer mixing ratio is determined based on the constraint that the global tracer mass be preserved (see Appendix A for details).The adjustment so obtained is small, with a mean value of −8 × 10 −5 ppm and a standard deviation of 0.015 ppm over a 1-year simulation.Since the surface pressure analyses are stored with only 16 bit precision (because the analysis uncertainty does not warrant more precision than this), this ultimately limits our knowledge of the tracer mixing ratio and the local tracer mass.With our global adjustment scheme, the adjustments are much smaller than 16 bit precision will allow, and thus the scheme is justifiable.

Model evaluation
We will use GEM-MACH-GHG to study the influence of uncertain atmospheric transport on the spatial scales recoverable in CO 2 fields, but, since this model has not yet been used for greenhouse gas transport, it is necessary to first document its ability to transport CO 2 .

Evaluation of meteorological fields
ECCC has been delivering operational weather forecasts for over 40 years and the products have been based on the GEM model for 19 years.As with many national weather forecast centres, ECCC participates in regular intercomparisons of operational forecasts following WMO standards.For the global configuration, forecasts on the medium (up to 10 days) range are compared every month and are available at http://web-cmoi.cmc.ec.gc.ca/verification/monthly/ observations/obs_monthly_e.html.Thus the quality of GEM transport has been documented.However, since changes were made to the model that affect weather prediction (specifically, the change to the continuity equation, the implementation of conservation of the dry air contribution to the global mean surface pressure and the convective transport scheme parameters) it is necessary to demonstrate the quality of the meteorological fields for the purpose of greenhouse gas transport.Since offline transport models often use reanalyses such as ERA-Interim or MERRA to transport constituents, we compare our 24 h forecasts against these products.In general, analyses more closely match observations than do forecasts because models are imperfect, and reanalyses should be superior in quality to operational analyses because the former use a superset of observations and a single, recent model version.Thus, our 24 h forecasts cannot be expected to be superior to any of ERAI, MERRA or JRA-55 (Kobayashi et al., 2015) reanalyses.Nevertheless, because of the usage of a 24 h forecast cycle (Fig. 1), it is important to verify that transport on this forecast range is reasonable, and reanalyses serve as high-quality reference fields.
The initial meteorological fields used at the start of every cycle (Fig. 1) come from archived operational products interpolated to our lower-resolution grid and topography.Thus our transport can only drift from reality for 24 h before it is corrected so that, even with our modifications to GEM, significant transport errors are not expected.Nevertheless, since changes of several hPa in 3-day forecasts of surface pressure do arise from the adjustment to the continuity equation and these could impact synoptic-scale forecasts, it is worth verifying the quality of our modified transport.Figure 2 (left column) shows the monthly mean difference of temperature, zonal and meridional wind differences between GEM-MACH-GHG and ERAI fields for July 2009 based on 6-hourly difference fields.ERAI fields were obtained from http://reanalysis.org on pressure levels in GRIB format and interpolated to GEM's grid at 1.5 • resolution.GEM fields were also output on the same pressure levels and resolution.The mean differences can be compared to those between other reanalyses: MERRA and ERAI (middle column) and JRA-55 and ERAI (right column).The latter difference fields were obtained from http://reanalysis.org using the Web-based Reanalysis Intercomparison Tools (WRIT; Smith et al., 2014) with the values saved in NetCDF format for plotting.It is evident that for all fields, GEM-MACH-GHG is as similar to ERAI as other reanalysis products although our products are 24 h forecasts (not analyses or reanalyses).For temperature (top row), the largest differences appear in the stratosphere for all models.In the stratosphere where models and observations are biased, assimilation can be challenging (Polavarapu and Pulido, 2016), so this result is not surprising.The extratropical zonal mean stratosphere is dominated by the slow Brewer-Dobson circulation, which is forced by waves propagating upward from the troposphere (Andrews et al., 1987;Vallis, 2006).The spectrum of waves from a given model depends on parameterizations such as convection and gravity wave drag which generate high-frequency subgrid-scale waves, as well as resolved waves.Therefore, different models can be expected to have rather different spectra and thus different forcing of the Brewer-Dobson circulation.For wind fields, the largest differences occur in the tropics.Zonal wind differences are greatest in the tropical stratosphere (middle row).The zonal mean tropical stratosphere is dominated by the quasi-biennial oscillation (QBO), which is driven by vertically propagating waves.Climate models have difficulty capturing the QBO as it depends on vertical resolution, gravity wave drag parameterizations and the spectrum of waves generated by tropical convection schemes (Baldwin et al., 2001;Campbell and Shepherd, 2005).Analyses can capture the QBO by assimilating the few radiosonde wind observations available in the tropics.However, the amplitude and phase of the QBO captured may depend on the underlying model and assimilation system characteristics.In general, since few direct measurements of winds are available to constrain analyses, and there are no simple dynamical balances available to infer winds from observations related to mass fields, model biases in the tropics are difficult to correct with data assimilation (Polavarapu and Pulido, 2016).In Fig. 2, the meridional winds (bottom row) differ most in the tropical troposphere where these issues of lack of measurements, and simple dynamical balances will prevail.
In summary, the meteorological fields used to transport CO 2 in our system are similar in quality to those from reanalyses which have been used with offline transport models.Reanalyses are not perfect and their uncertainty can be quantified through the degree of disagreement with other reanalysis products.Thus, when offline transport models use reanalyses, this uncertainty impacts their transport quality.While we use a 24 h update cycle, this subsection has demonstrated that throughout the 24 h forecast the uncertainty is comparable to that of an individual reanalysis product.Thus results on the sensitivity of transport error to meteorological uncertainty obtained with our sequence of 24 h forecasts has direct relevance to transport errors obtained with a 6 h insertion of reanalysis products.While statistics for July 2009 were shown in Fig. 2, those for December 2009, July 2010 and December 2010 are found in the Supplement (Figs.S3-S5).

Evaluation of CO 2 fields
Having established that the meteorological fields that will transport constituents are sufficiently accurate, the CO 2 evolution can be assessed.Since the goal is to assess GEM-MACH-GHG as a transport model, realistic surface-toatmosphere fluxes are required.Without realistic sources and sinks of CO 2 , discrepancies with observations could be equally attributed to erroneous fluxes as to erroneous transport.Even with a good source of retrieved atmospheric fluxes, the evaluation of GEM-MACH-GHG transport requires care (as discussed below).For this purpose, posterior fluxes from CarbonTracker 2013B (hereafter referred to as CT2013B; Peters et al., 2007) were obtained and regridded in a mass conservative way to GEM's grid.CarbonTracker was chosen because it is generally recognized as a good product, is regularly monitored and updated and is readily available from http://carbontracker.noaa.gov.These fluxes were also used for a similar purpose in Houweling et al. (2010).CT2013B fluxes are available every 3 h, with fluxes inserted at every model timestep.With an initial condition from Car-bonTracker for 1 January 2009 00:00 UTC and the CT2013B fluxes, GEM-MACH-GHG was run for all of 2009-2010.The resulting CO 2 simulations are compared against observations assimilated by CT2013B (continuous observations), observations not assimilated by CT2013B (columns from Total Carbon Column Observing Network (TCCON) and aircraft profiles) and gridded 3-D concentration fields from CT2013B. Figure 3 compares CO 2 time series from GEM-MACH-GHG with CT2013B fluxes to surface observations from ECCC's greenhouse gas measurement network (Worthy et al., 2005) at Alert, East Trout Lake and Sable Island.Alert is a remote Arctic site far from CO 2 sources, so it can be used to assess long-range transport.East Trout Lake is close to sources and will reflect diurnal variations in CO 2 fluxes The mean diurnal cycle at individual stations can reveal more clearly the realism of a model's boundary layer variation.Here we choose two sites close to sources and sinks because the amplitude of the diurnal cycle at such sites should vary through the year.Figure 4 shows that the model's mean diurnal cycle at East Trout Lake with GEM-MACH-GHG with CT2013B (blue) fluxes compare well to measurements (black), whereas CT2013B (green) has a too-large amplitude in April and July 2009 (consistent with Fig. 3).However, the model's behaviour varies with location and time.At Fraserdale, GEM-MACH-GHG's diurnal cycle amplitude is clearly too low in all months whereas CarbonTracker fares better in July and October 2009.The fact that the diurnal variability from GEM-MACH-GHG is lower than that observed stems from the choice made for the minimum value of eddy diffusivity (Sect.3.5).Lowering this value can increase the amplitude of the diurnal cycle, but with the increased risk of occasional spuriously low values of CO 2 during summer daytime.As noted earlier, it is difficult to compare models to measurements on sub-diurnal timescales (Law et al., 2008;Patra et al., 2008) and most models have difficulty in capturing boundary layer evolution so flux inversions typically use only afternoon mean measurements.
GEM-MACH-GHG transport is directly compared to Car-bonTracker's transport in Figs. 5 and 6. Figure 5 presents the column mean CO 2 weighted by air mass on 1 July 2009 00:00 UTC and 31 December 2009 00:00 UTC from GEM-MACH-GHG with CT2013B fluxes (top panel), CT2013B (middle panel) and the difference between these two (bottom panel).The higher resolution of GEM-MACH-GHG is evident particularly in July, when the largest differences occur at small spatial scales.CT2013B has more CO 2 than GEM-MACH-GHG in both the tropics and the Northern Hemisphere.However, in winter (right column) CT2013B has more CO 2 mainly in the tropics.The tropical differences in both seasons are likely related to differences in convection schemes in the two models.The increased CO 2 in the northern midlatitudes in boreal summer may be due to differing meridional transport or to differing rates of vertical mixing or a combination of both.Zonal mean fields also reveal that the greatest differences between the two GEM-MACH-GHG simulations and CT2013B are in the summer (Fig. 6, top panel).Throughout the troposphere and near the surface, CT2013B has more CO 2 .In winter, CT2013B has a slight deficit near the surface compared with GEM-MACH-GHG.Differences in the stratosphere between the two models are also evident in Fig. 6.Since GEM-MACH-GHG has better vertical resolution compared to CarbonTracker (80 vs. 25 levels) and GEM is designed to have a realistic stratosphere (Charron et al., 2012), differences in stratospheric and mesospheric flow are to be expected.However, the mass of CO 2 in the stratosphere and mesosphere is very small, so the column mean or surface values would be insensitive to such differences.In summary, when GEM-MACH-GHG is run with CT2013B fluxes, differences in transport errors between CarbonTracker and GEM-MACH-GHG are evident.Differences are within 3 ppm in column mean and 4 ppm in zonal mean.Neither (any) model can be expected to have perfect transport, so the acceptability of transport is generally gauged through comparisons of model predictions to measurements.
To assess seasonal timescales, it is useful to compare to the TCCON measurements (Wunch et al., 2011).The data used for this study are from the GGG2014 release, available on the network's website http://tccon-wiki.caltech.edu.All sites with measurements in 2009 and 2010 are selected, as listed in Table 1.The details of how the model profiles were converted to column-averaged dry mole fraction (X CO 2 ) and smoothed following Wunch et al. ( 2010) is provided in Appendix B along with the precise definitions of the statistics discussed here (bias, root mean square (RMS) and scatter).The statistics for 2009-2010 are shown in Table 2.The  bias is below 1 ppm for every station except Eureka, which only had 49 h of measurements in 2010.No selection is applied when considering which sites are included in the "all", "mean" or standard deviation (SD) statistics.At Eureka, the −2.66 ppm bias significantly impacts the station-to-station SD of the bias, it is 0.3 ppm without Eureka.(Note that an error in surface pressure was recently discovered at Eureka and preliminary results for a correction suggest a 0.5 ppm difference in X CO 2 , which would reduce the Eureka EC-CAS bias by a corresponding amount.)Except for Eureka and Park Falls, the station bias is positive and the overall standard deviation is 0.96 ppm with a high correlation coefficient of 0.95.Although we did not perform a data assimilation, the posterior fluxes from CarbonTracker contain information from the observations they used, so we can compare our Table 2 to Massart et al. (2016, Table 2).The biases of individual stations are mostly lower here, as is the overall averaged bias.This is because the surface observations assimilated by CarbonTracker with a long assimilation window are able to constrain the global atmospheric growth rate, whereas the system used by Massart et al. (2016) does not use long assimilation windows and thus does not constrain Table 2. Statistics for the average hourly X CO 2 comparison between TCCON measurements and GEM-MACH-GHG simulations: RMS (ppm), bias (ppm), scatter (ppm) and correlation coefficient R. The "all" line shows the mean of each parameter using data from all sites combined."Mean" is the average of each parameter and SD is the standard deviation of the station's bias.N is the number of data pairs (or sites for mean and SD) used in the computation of the statistics.the global growth rate as effectively.The full time series for Park Falls (USA) and Wollongong (Australia) are shown in Fig. 7 and the seasonal bias as well as the seasonal statistics using data from all sites combined are shown in Table 3. Seasons are defined as DJF (December, January and February), MAM (March, April and May), JJA (June, July and August) and SON (September, October and November).The model is able to reproduce seasonal variations of X CO 2 with biases ranging between −1 and +1 ppm (excluding Eureka) and scatter values consistently below 1 ppm.The large negative bias at Eureka in autumn is consistent with the time series shown in Fig. 3. Other northern stations (Sodankylä and Bialystok) have similar but smaller biases.As discussed earlier, the discrepancy of GEM-MACH-GHG transport with that of CarbonTracker (as imprinted in the posterior fluxes) may explain this behaviour since other posterior fluxes do not have this particular issue (Fig. S6).

Site
The vertical structure of model CO 2 is compared to NOAA aircraft profiles (Sweeney et al., 2015) over Canada and the USA in Fig. 8.Following Agustí-Panareda et al. (2014), mean model profiles at the nearest model grid point and timestep to the profile location were averaged over all profiles for a season.Both observed and model values were binned into 1 km layers.The observations are from Ob-sPack2013 (2013) (Masarie et al., 2014) and include only profiles from Canada or the contiguous USA.The annually averaged model profiles are shown in panel a.In the annual average, GEM-MACH-GHG has good agreement with these independent measurements while CT2013B has a very slight positive bias.Both models are quite good compared to the ensemble of models shown in Stephens et al. (2007, Fig. 2b).However, the other panels reveal that GEM-MACH-GHG's excellent annual result in the free troposphere is because of compensating errors in different seasons.The boreal winter season (DJF) is not shown because December 2008 was not simulated.However, the behaviour of the model in boreal winter 2009 (based on January-February) and 2010 is qualitatively similar to its behaviour in spring (panel b).Panels b-d reveal that GEM-MACH-GHG agrees quite well with observations from 3 to 6 km in all seasons.However, from In autumn (panel d) the gradient is slightly too large, while CT2013B is almost perfect.Thus, the vertical mixing just above the boundary layer is too weak, as with most models (see Yang et al., 2007).A possible reason for the overestimation of vertical gradient with GEM-MACH-GHG is the lack of tracer transport through shallow convection.These biases in vertical gradients will be relevant for regional flux inversions (see Stephens et al., 2007) that may use GEM-MACH-GHG results as boundary conditions.atmospheric analyses Having established that GEM-MACH-GHG can simulate CO 2 reasonably well, we turn our attention to the question of how atmospheric transport modulates the CO 2 distribution.The evolution of CO 2 can be described by the species transport equation and thus may be considered to be perfectly predictable.However, this is only true if the advecting fields are perfectly known, and this is never the case.With a coupled meteorology and forecast model, the impact of the uncertainty of meteorological fields on CO 2 transport can be explored.In a data assimilation or flux inversion system, when the fluxes are well constrained by observations, the ability to estimate the CO 2 fields using observations will ultimately be limited by the loss of meteorological predictability, just as the quality of weather forecasting products are.Thus it is useful to identify these limits on the spatial scales that can be retrieved in analyses of CO 2 , because these predictability lim-its can then be compared to the spatial scales of CO 2 retrievable in the presence of meteorological analysis uncertainty.In this section, we use GEM-MACH-GHG to determine the predictability of CO 2 on weather and climate timescales with the latter referring to subseasonal to seasonal scales.In other words, we isolate and study the component of transport error due to meteorological state (initial condition or analysis) errors.In Sect.5.1, the classic weather predictability problem of uncertainty in meteorological initial conditions is considered, whereas longer timescales and the transport errors due to uncertain meteorological analysis errors and model errors are considered in Sect.5.2.

Weather timescales
Although flux inversion systems focus on retrieving relatively long timescale signals (than 2 weeks), it is useful to first consider CO 2 predictability on weather timescales before considering errors o longer timescales (next subsection).Specifically, we isolate the impact of the loss of meteorological predictability on CO 2 predictability.The meteorological predictability problem on weather timescales is related to forecast sensitivity to initial conditions, but the atmospheric variability of CO 2 on diurnal (Law et al., 2008), synoptic (Chan et al., 2004;Agustí-Panareda et al., 2014) and seasonal and interannual (Gurney et al., 2002(Gurney et al., , 2004;;Baker et al., 2006a;Le Quéré et al., 2015) timescales is largely governed by the terrestrial biospheric fluxes and hence is determined by sensitivity to boundary rather than initial conditions.Nevertheless, the CO 2 predictability error arising from loss of meteorological predictability on weather timescales has not to our knowledge been identified, and it can be used to iden- tify an upper limit on forecast errors.This may be relevant for operational data assimilation or forecasting systems such as those at ECMWF (Agustí-Panareda et al., 2014) and NASA Goddard (Ott et al., 2015), which use update cycles of 12 or 24 h and also examine the quality of short-term forecasts.It will also serve as an upper limit for transport error arising from the presence of uncertain meteorological analyses in Sect.5.2.
Predictability of weather normally refers to the sensitivity of forecast errors to initial conditions such that any infinitesimal perturbation will lead to diverging forecasts in a finite length of time.This is the so-called butterfly effect and it occurs because of the underlying nonlinear chaotic dynamics of the governing equations (Palmer, 2006).To compute predictability error of meteorological variables on weather timescales, one can simply start with a reference simulation and perturb the initial conditions.Eventually, the forecasts will diverge, but the error will saturate at climatologi-cal levels.Once saturation has been reached, the statistics of this predictability error can be determined.However, with a transported tracer such as CO 2 , the model forecast requires regular insertion of wind fields.During our forecast cycle, the meteorology is constrained by analyses every 24 h and departure from reality will represent at most a 24 h forecast error.Thus for transported constituents, the definition of the predictability experiment is slightly different.The reference simulation will be taken as the GEM-MACH-GHG 2-year run with CT2013B fluxes.Then a comparable "climate cycle" is run in which the model, initial conditions, CO 2 fluxes and surface forcing are identical to those used in the reference cycle.However, with the second and all subsequent cycles, the meteorological fields are not replaced by analyses but are instead copied from the 24 h forecast fields.Thus the meteorology fields used to transport the CO 2 field in this "climate cycle" are never updated with observations (analyses) and will thus depart from those used in the reference cycle Table 3. Seasonal bias (ppm) for the hourly averaged X CO 2 comparison between TCCON measurements and GEM-MACH-GHG simulations.The mean (using data from all sites) bias (ppm), scatter (ppm) and correlation coefficient (R) are also shown under "all".Seasons with fewer than 10 pairs are not included in the "all" calculations.

Site
Latitude in the first 2 weeks.The divergence of the CO 2 field in the "climate cycle" from that in the reference cycle, once the error has saturated, defines the CO 2 predictability error arising from the loss of meteorological predictability.
Starting from 1 January 2009 00:00 UTC, the reference and climate cycles are run for 1 month with fields saved every 6 h.The differences between the corresponding fields from the two cycles are computed for temperature, zonal and meridional wind components, CO 2 and surface pressure.Vorticity and divergence fields are computed from the wind difference fields.Then the global mean of the zonal standard deviation of each difference field is computed after first subtracting the zonal mean.The resulting global mean values are called predictability errors and they have the same units as the corresponding forecast field.In order to get a comparable scale for all variables, the errors were normalized by the global mean of the zonal standard deviation for a reference state.The choice of this state is arbitrary but has implications on the maximum values attained (as will be discussed below).The reference field was taken as the initial state used to launch both the control and climate cycles and corresponds to 1 January 2009 00:00 UTC. Figure 9 shows the time and height variation of the predictability error normalized by the variability of the reference state for all four variables.When the predictability error approaches the variability of the reference state, values approach 1.Thus predictability is expected only when the relative error is much less than 1.In Fig. 9, we see that temperature loses predictability within 10 days, as expected (i.e. the normalized error reaches 0.8).However, CO 2 loses predictability in the troposphere within 2-3 days except very near the surface, where it reaches 5 days.The predictability of CO 2 more closely resembles that of the wind (vorticity and divergence) fields which also lose predictability in fewer than 5 days in the troposphere.This makes sense because the wind fields are used to transport the CO 2 field.The difference in evolution of the CO 2 field in the reference and climate cycles (figures not shown) reveals largest values to be associated with gradients created by large fluxes (whether natural or anthropogenic).This ability of the uncertainty in wind analyses to act on CO 2 gradients and to spread the uncertainty downstream was previously illustrated by Liu et al. (2011) using an ensemble of wind fields.CO 2 is more predictable in the stratosphere, with the loss of predictability occurring after 5 days in the lower stratosphere.The extended predictability of CO 2 in the lower stratosphere is similar to that seen in the vorticity field (bottom left panel).The reason that the vorticity field is more predictable than the divergence field (with a loss of predictability occurring after 3-4 days in the troposphere) is because the vorticity field is associated with slower rotational modes whereas the divergence field is often associated with higher-frequency waves.The atmospheric kinetic energy spectrum is dominated by rotational motions in the troposphere (Koshyk et al., 1999;Skamarock et al., 2014).In the stratosphere, the zonal mean flow in winter is driven by very large-scale vertically propagating planetary waves (Andrews et al., 1987;Vallis, 2006) so large-scale rotational modes dominate the energy spectrum and extended predictability in vorticity and CO 2 results.While predictability is considered lost when its normalized error approaches 1, sometimes the relative error is much greater than 1.This occurs because of the arbitrary choice taken for the normalization.In our case, the reference state corresponds to the initial state, which corresponds to a relatively quiescent synoptic situation since January 2009 marks the strongest and most prolonged stratospheric major warming on record (Manney et al., 2009).The criteria for a stratospheric sudden warming (SSW) were met on 24 January 2009 when zonal mean easterlies replaced the climatologically normal westerlies at 10 hPa.However, easterlies were noted in the mesosphere prior to this date (Manney et al., 2009).This pattern is consistent with the appearance of anomalously large predictability error in relative vorticity in the mesosphere prior to 24 January and the appearance of anomalously large errors in the stratosphere (layer 6) after 15 January.Although global mean values are shown, it is zonal standard deviation that is computed, and the departure from a zonal mean will be large during a wave 2 vortex splitting event such as occurred in 2009.Thus zonal standard deviations are anomalously large throughout much of the Northern Hemisphere because the climate cycle does not capture this event whereas the reference cycle does.The extent of the disturbance in the Northern Hemisphere is large enough to influence the global mean values.Anomalously large CO 2 predictability error appears in the mid-stratosphere around 20 January 2009.The disturbance of CO 2 due to the SSW is to be expected since the disturbance of CO in the mesosphere, N 2 O in the mid-stratosphere and H 2 O in the lower stratosphere as the vortex deformed and split was evident in MLS (Microwave Limb Sounder) observations (Manney et al., 2009).
Figure 10 compares the layer mean normalized predictability error for different variables for the layers indicated by the dashed lines in Fig. 9 and labelled in white text at the left edge of each panel.Near the surface, the normalized CO 2 predictability error closely follows that of the specific humidity field (for about 3 days).Both moisture and CO 2 fields are advected by the wind fields and are similarly affected by the predictability of the wind fields.In this bottom layer, CO 2 and moisture are more predictable than the wind field presumably because of their dependence on surface fluxes, although the normalization, which is layer dependent, cannot be ruled out as a contributing factor.However, they are both less predictable than the temperature field (Fig. 10a).In the troposphere, the loss of CO 2 predictability is similar to the loss of predictability in vorticity (panels b, c) for the first few days.After that, the loss of predictability for CO 2 is faster than that for the vorticity and divergence fields.While the normalization may play a role in this result, it is also worth noting that the predictability of vorticity and divergence fields will differ from that of wind components because very smooth wind field errors generate little error in vorticity and divergence fields.Throughout the atmosphere, temperature loses predictability at a slower rate than CO 2 does.While both moisture and CO 2 are transported by wind fields, moisture is also a dynamic variable; thus the loss of predictability for specific humidity is not the same as that for CO 2 in the lower to mid-troposphere (panels b, c) beyond day 5. Predictability error increases from day 1 to reach saturation levels in 10-15 days for all levels except the upper stratosphere (Fig. 10f), which is affected by the SSW.Thus in the next section, climatological levels of predictability error will be discussed for periods longer than 1 month in the troposphere.
Waves of 24 h in period are seen, particularly, in the vorticity and divergence plots in Fig. 10.This occurs because the forecasts in the reference cycle are abruptly corrected every 24 h with the insertion of a new analysis.When the predictability error of the 24 h forecast error is large compared to that of the analysis valid at the same time, we see stripes at the 24 h period.Thus we conclude that the normalized 24 h forecast error of the wind field is much larger than that of the temperature field.This makes sense because the global observing system is dominated by information about the mass field with relatively sparse direct observations of the wind field (Baker et al., 2014).In addition, the mass field (which is reflected in the temperature field) is a much smoother field and is thus more easily observable with a given network rel- . Time series of layer mean normalized predictability error during January 2009 for CO 2 (black curves), temperature (red curves), vorticity (blue curves), divergence (green curves) and specific humidity (cyan curves).For clarity, no curves for specific humidity are plotted in the bottom row of panels.The normalized errors are averaged over 12 model levels; the top and bottom levels used in the average are given in approximate pressure above each panel.The layer numbers are associated with the layers defined by dashed horizontal lines in Fig. 9 and correspond to the near surface (layer 1, a), the lower troposphere (layer 2, b), the mid-troposphere (layer 3, c), the upper troposphere (layer 4, d), the lower stratosphere (layer 5, e) and the upper stratosphere (layer 6, f).
ative to fields which are dominated by smaller spatial scales (such as vorticity or divergence).
In summary, the global predictability of CO 2 due to uncertain meteorological initial conditions is very short in the free troposphere and is associated with the predictability of wind fields.This predictability limit refers only to sensitivity to meteorological initial conditions and it will be counter balanced by the predictability coming from biospheric fluxes on diurnal and synoptic scales.It is important to note (and this is discussed more fully in Sect.6) that the predictability diagnostic used here is dependent on an arbitrary normalization and an arbitrary threshold to define predictability so absolute results are not expected.Instead, the relative predictability between variables or between atmospheric layers is expected to be more reliable.To improve this type of predictability, more observations will be needed where the wind fields have finer spatial scales and where convection is occurring.The current global meteorological measurement network is relatively sparse in the tropics where convection is important, but new observations from space-borne lidars may be able to remedy this problem (Baker et al., 2014).

Seasonal timescales
As noted earlier, CO 2 observations contain information on seasonal to interannual timescales.Specifically, the global surface network used in flux inversions is able to constrain the global CO 2 budget and capture seasonal and interannual variability of the global fluxes (e.g.Baker et al., 2006a;Peylin et al., 2013).The source of predictability on subseasonal to seasonal and longer timescales partially derives from climate predictability on those scales but also from long timescale information contained in terrestrial biospheric fluxes which are, in turn, influenced by climate variability (Patra et al., 2005).In this section, we explore the predictability of CO 2 on longer timescales and compare these to CO 2 simulation errors due to the use of uncertain meteorological analyses.
The predictability experiment in Sect.5.1 represents an extreme case in which no information from observations is present in the wind fields after the initial time.In reality, in our CO 2 transport cycle (Fig. 1), the wind fields are constrained to observations by the insertion of a meteorological analysis every 24 h.However, the analyses are not perfect and have a certain level of uncertainty.To simulate this uncertainty, we could perturb the analyses every 24 h with an analysis error.In a variational data assimilation system such as that used at ECCC, it is possible to estimate the analysis error covariance matrix but it is expensive to do so, and such estimates are not routinely made.In contrast, a simple perturbation such as random spatially uncorrelated errors will not be useful as they will primarily generate unbalanced motions.What is more relevant is a perturbation of the size and shape of the 6 h analysis error, given the use of a 6 h forecast cycle in operations.(Reanalyses are also available at 6 h intervals and these are sometimes used to constrain flux inversions.)Thus, in order to simulate a coherent 6 h analysis The lower and upper model levels averaged are indicated above each frame.Approximate pressure is obtained from model level by multiplying by 1000 (which corresponds to assuming a reference surface pressure of 1000 hPa).The CO 2 reference state spectra (blue curves), predictability error (black curves) and error due to a 6 h shift in analysis fields (red curves) are shown.
error, we simply insert the analysis state valid 6 h prior to the actual analysis time (i.e. the one from 18:00 UTC of the day before, instead of the correct one from 00:00 UTC), relabelling the date and time to the correct ones.In addition, because the shift in the diurnal cycle for the meteorology would impact the CO 2 predictability, the diurnal cycle is removed from the perturbations by subtracting the monthly mean of each synoptic hour.Then, the deviation of the CO 2 field from this perturbed analysis cycle from the reference cycle defines the error due to the use of uncertain meteorological analyses.This error should be much smaller than the predictability error arising from uncertain meteorological initial conditions.However, it should be larger than an actual analysis error because of the additional component corresponding to the evolution of the true state in 6 h.Thus it is an overestimate of analysis error.
Figure 11 shows the monthly mean spatial spectra of various difference fields averaged over several model levels for July (top row) and December (bottom row) of 2009.The spectra refer to the spherical harmonics (Boer, 1983) of a scalar field multiplied by its complex conjugate and summed over zonal wavenumbers.The x axis then defines a total wavenumber.While spectra were computed for each day at 00:00 UTC, these were averaged over the month to filter some noise and identify a robust signal.In addition, they were averaged over 12 model levels to get representative spectra for a few atmospheric layers, namely the bottom 4 layers shown in Fig. 9.The blue curves in Fig. 11 depict mean spectra of the CO 2 state from the reference cycle similarly averaged in time and in the vertical dimension.The black curves represent the predictability error arising from uncertain meteorological initial conditions.This error is very small at the end of the first cycle (on 2 January 2009) but rapidly increases during the first 2 weeks to saturate at climatological levels (Sect.5.1).Since we are interested in this saturated level of error, we do not consider the first month of errors.The 2 months chosen in Fig. 11 represent the variation seen in various months of the year.The predictability error is seen to be lower than the reference state itself for very large spatial scales but quickly equals (around wavenumber 10) and then surpasses the power in the reference state.The reason that the power in the predictability error can be larger than that in the state itself is that it involves the difference of two fields.In the limit where two fields become uncorrelated, the variance of the difference equals the sum of the variances.If the two fields have the same climatological variance, the variance of the difference is twice the climatological variance.Thus it is not surprising that the power in the predictability error should surpass that in the reference state for small spatial scales.What is more intriguing is that some information is still retained in the largest scales (wavenumbers less than 5) even after 6 or 12 months of simulation.The source of this predictability at very large scales is mainly due to surface fluxes of CO 2 which are the same for both the reference and climate cycles.However, it might also be partially from the surface forcing of meteorological fields.Subseasonal to seasonal predictability is manifested in modes of variability such as the Madden-Julian oscillation (MJO), the Pacific North American (PNA) pattern, midlatitude blocking events and the North Atlantic Oscillation (NAO) (Waliser, 2006) and their predictability derives from atmospheric boundary conditions, namely sea surface temperature, soil moisture, snow cover, vegetation and sea ice (Shukla and Kinter III, 2006).These ocean and land surface conditions influence fluxes of moisture and sensible and latent heat into the atmosphere which may change low-level atmospheric convergence and lead to atmospheric heating anomalies, which influence the large-scale flow.To see if long timescales in the meteorological analyses play a role in the large-scale predictability seen in Fig. 11, an experiment was run in which the predictability experiment was repeated, but this time incorrect surface fields (from 3 months later) were used.With no information from atmospheric observations as well as a seasonally shifted error in the surface forcing, the predictability error is worsened at these largest scales, particularly in the summer (June, July and August) near the surface (Fig. S7) and in the lower troposphere (not shown).The differences in CO 2 evolution in the two predictability experiments during boreal summer are largest in the northern extratropics (not shown).This confirms that the ocean and land surface are playing a role in predictability of the system at the largest scales in the lower troposphere in boreal summer.In the mid-and upper troposphere, smaller impacts are seen but the impact is largest in the spring (not shown).In Fig. S7, the remaining predictability seen at wavenumbers below 10 for all months is then attributed to the common CO 2 fluxes used by the reference and predictability experiments.Thus Fig. S7 confirms that the dominant source of predictability at large scales seen in Fig. 11 is due to surface fluxes.
In summary, a direct impact of climate predictability on CO 2 predictions through the ocean and land surface is seen through the worsened predictability in boreal summer months when biospheric CO 2 fluxes are largest.Since CO 2 fluxes were specified, and were the same in the reference and predictability experiments, this climate signal is retained and explains most of the CO 2 predictability at large scales in Figs.11 and S7.Finally, it is worth noting that a shift of 1 month in surface fields resulted in no real deterioration of the predictability error since there is still significant correlation among surface fields due to a 1-month lag.
The red curves in Fig. 11 depict the spectra due to analysis uncertainty.As expected, these spectra are reduced compared to predictability errors because meteorological analyses are used in this cycle, although they are 6 h out of date.In particular, significant error reduction at large scales is seen.However, the red curves also intersect the reference state spectra at increasingly smaller wavenumbers as height increases.Thus, if analysis uncertainty is considered, there is a gain of information over climatological error levels defined by the predictability error but only for the larger spatial scales.The analysis perturbations used here overestimate analysis error so we can expect that the spatial scales resolved are a conservative estimate.However, whatever the size of the error, there will be a limit to the spatial scales resolved because of meteorological analysis imperfections.At the very least, the meteorological observing system has a given spatial resolution and, moreover, meteorological assimilation systems (like ECCC's) may compute analysis increments at a coarser resolution than that of the model.Beyond the point where the red and blue curves intersect, there is no useful information in the CO 2 field at these scales due to meteorological analysis uncertainty because the power in the CO 2 prediction error is larger than that in the CO 2 state.In fact, the CO 2 predictions with analysis uncertainty (red curves) asymptote to the predictability error spectra for large wavenumbers.For these spatial scales, using the updated meteorological analysis is no better than having no updates at all.Near the surface, there is the greatest gain of information, but in the upper troposphere the spectra of CO 2 errors due to analysis uncertainty are less than that of the reference state only for wavenumbers lower than 30.Thus the fact that the meteorological analysis has information on only certain spatial scales places limits on the spatial scales that can be retrieved in a transported field such as CO 2 .
Another view of the spectra can be obtained by summing over total wavenumber and plotting with respect to zonal wavenumber (Fig. 12) because zonal wavenumber spectra are more indicative of the tropical signal.In this figure, it is evident that there is almost no information retained in the predictability error (i.e. the black curve lies above the blue one).In other words, predictability error spectra surpass that of the reference state at all levels except near the surface for the first few zonal wavenumbers.When a 3-month shift in surface fields is included in the predictability experiment, all predictability is lost in July and August since the power in the predictability error (red curves) exceeds that of the reference state (black curves) for all wavenumbers (Fig. S8).Predictability is mostly lost in September as well.Thus using the correct land and ocean surface fields may be relevant for capturing predictability in July, August and September in the tropics.The dominant mode of tropical subseasonal variability is the MJO, which is characterized by very large scales (zonal wavenumbers 1 and 2 in wind and rainfall fields) in the tropics (Waliser, 2006); GEM can capture predictability associated with this mode (Lin et al., 2008), supporting the notion that GEM has predictive skill in the tropics on seasonal scales.However, since predictability was already close to lost even with the correct land and ocean surface (power in predictability error is same order of magnitude as that for reference state), the worsened predictability in July to September with incorrect land and ocean surface fields is not a clear indication of their influence in the tropics.Indeed, outside of these 3 months, shifting the land and ocean surface fields by 3 months has little impact (red and black curves are similar in Fig. S8).This suggests that most of the information retained at the largest scales (seen in Fig. 11) in the predictability experiment is coming from the extratropics or signals with latitudinal structure for these months.Indeed, the CO 2 evolution (not shown) in the reference and predictability experiments differs the most in northern hemispheric extratropics.Thus the large-scale predictability seen in the northern extratropics is due to large CO 2 variability during October to May associated with biospheric fluxes, but from June to September the ocean and land forcing of GEM's climate is also important.When analysis uncertainty is simulated (red curves), there is a gain of information over predictability error but fewer than 20 zonal waves are resolved in the mid-troposphere (Fig. 12, layers 2 and 3).In the lowest layer, about 40 waves are resolved.Compared to Fig. 11, there are significantly fewer waves being resolved.Thus structure in the meridional direction is better resolved in analyses than structure in the zonal direction is.
It is also useful to examine other sorts of model errors in terms of their impact on CO 2 predictions.Figure 13 shows the spectra of errors due to the inclusion (or not) of convective tracer transport (cyan curves).The impact of adding convective tracer transport is primarily at large scales and exceeds CO 2 errors due to imperfect wind analyses (red curves) only in the mid-troposphere for wavenumbers less than 5. Spectra for April 2009 were shown in Fig. 13 because the impact of adding convective tracer transport was largest in spring months.The impact is always less than that due to imperfect wind analyses if zonal wavenumber spectra are considered (bottom row).Thus, for our model, this type of model error is exceeded by transport error due to uncertain meteorological analyses for most spatial scales.This type of analysis in which components of transport error are compared may thus be useful in diagnosing the dominant sources of transport error for a given transport model.

Summary and discussion
A new capability for simulating CO 2 using ECCC's operational weather and environmental prediction models has been developed.The adaptations required for greenhouse gas simulation include the implementation of a global mass fixer for the semi-Lagrangian tracer transport scheme, the implementation of a mixing ratio defined with respect to dry air for tracer variables, the addition of convective tracer transport and modification of a parameter in the boundary layer scheme.A sequence of 24 h meteorological forecasts is used to transport CO 2 fields in a forecast cycle involving a coupled meteorological and tracer transport model.The 24 h meteorological forecasts are as similar to the ERA Interim reanalyses as other reanalyses are .That means that throughout the sequence of 24 h forecasts, the meteorological uncertainty is comparable to that of a reanalysis dataset which could be used to constrain an offline transport model.Using prescribed posterior fluxes from NOAA's CarbonTracker (CT2013B), the transport of the model has been assessed.The CO 2 fields compare well to observations assimilated in the posterior fluxes (surface hourly measurements) as well as to independent observations (within 2 ppm for TCCON and NOAA aircraft profiles) and to CarbonTracker mole fractions (within 3 ppm in column means).Synoptic and seasonal timescales are well captured but, as with most transport models, the diurnal cycle ampli-tude is too low in summer.The vertical gradient in the midtroposphere is slightly overestimated but the gradient from 1 to 3 km does not agree as well with observations and the error in the gradient changes with season.
The predictability of CO 2 forecast uncertainty.CO 2 forecasts or are implicitly used in flux inversions when computing the model-data mismatch.The forecast or prediction error of a transport model is often referred to as "transport error" and transport error is a significant component of posterior flux errors from flux inversions.Transport error is comprised of meteorological state, CO 2 state, flux and model formulation errors.While considerable research has been devoted to better understanding the nature of model formulation and flux errors on CO 2 predictability (or CO 2 transport error), much less is known about the nature of meteorological state errors and their impact on CO 2 transport errors.Liu et al. (2011) have shown that transport errors due to meteorological uncertainty are 1.2-3.5 ppm near the surface and 0.8-1.8ppm in the column mean.Here we examine the atmospheric processes and scales implicated in such impacts.Of course, flux and model formulation errors are a significant component of transport error and the various contributions to transport error do not act in isolation.However, isolating this impact provides new insight into transport error processes.Moreover, it permits the identification of the spatial scales of CO 2 that can be trusted in the context of meteorological analysis uncertainties.
The fact that weather becomes unpredictable beyond 2 weeks has implications for CO 2 transport.With an experimental design developed to isolate the impact of weather predictability due to imperfect initial conditions on CO 2 , we have demonstrated that CO 2 also has limited predictability.That is, uncertainty in meteorological initial conditions is evident in CO 2 simulations when other sources of transport error are removed.This is a newly identified source of uncertainty which is primarily of pedagogical interest because transport models always have updated information from meteorological analyses.However, isolating this error process leads to insights into the processes which couple CO 2 and meteorological errors.Moreover, this type of error forms the upper limit of transport errors due to uncertain meteorological analyses.The predictability of CO 2 due to uncertain meteorological initial conditions is found to be short (just a few days) and is comparable to that of the wind field during this time.After a few days, CO 2 loses predictability faster than the wind field.Predictability is greatest near the surface and in the lower stratosphere.Reduced predictability in the stratosphere was seen during the prolonged and strong SSW of January 2009.Predictability of CO 2 is shorter than that for temperature at all levels but is similar to that of the specific humidity field for the first few days, after which CO 2 becomes less predictable than moisture.The predictability measure used is based on zonal variability.In addition, the measure is normalized by the variability of the initial state so that loss of predictability by different variables could be compared.Since there is an arbitrariness in the choice of the normalization state, and in choosing a threshold to define predictable states, absolute times for predictability cannot be obtained.Rather, it is the relative predictability that is expected to be more reliable.Specifically, the greater loss of predictability for CO 2 than the temperature, humidity or wind fields (the latter after a few days) and the increased predictability of CO 2 near the surface and lower stratosphere compared to that in the free troposphere are the new results.While these results are model specific (as are predictability results in general) and metric specific, there are reasons to expect more generality of the results.Increased predictability of CO 2 near the surface can be explained by the constraint of fluxes near the bottom boundary.Increased predictability of CO 2 in the lower stratosphere makes sense because of dominance of slow rotation motions in that layer of the atmosphere.Nevertheless, these results need to be confirmed by other models and metrics.Massart et al. (2016) have computed the predictability of X CO 2 forecasts obtained from analyses which assimilated GOSAT observations and found predictability up to 5 days globally.Results were neither produced for CO 2 at various altitudes nor compared with those for meteorological variables.However, the predictability of X CO 2 was computed for different regions.Predictability was greater in the tropics than in the extratropics.The measure used was an anomaly correlation wherein forecasts anomalies with respect to a climatology are correlated with analysis anomalies.This measure of predictability is useful for identifying predictability of a given forecast system at specific dates and forecast ranges.However, as with all metrics there are known deficiencies.In particular, the anomaly correlation can be optimistic where observation density is poor because forecasts and analysis anomalies will be perfectly correlated in the extreme case of no observations.Thus, the predictability results of Massart et al. (2016) are not directly comparable with those obtained here.
Even with no information from meteorological analyses, the predictability of CO 2 was seen on seasonal timescales for very large spatial scales and was found to be primarily due to surface fluxes.However, worsened predictability of CO 2 was seen on these scales in the northern extratropics in boreal summer if the surface atmospheric forcings were incorrectly specified.In other words, the land and ocean surface also play a role in CO 2 predictability on seasonal timescales.The proposed mechanism is that the land and ocean surface perturb atmospheric circulations which transport CO 2 and the impact is greatest where CO 2 gradients are largest: in the northern midlatitudes in boreal summer.
Our predictability error is defined by the error due only to imperfect meteorological initial conditions.Thus compared to a reference cycle, the perturbed cycle has no updates of meteorology with analyses.The more realistic situation of imperfect meteorological information was addressed by perturbing the meteorological analyses used in a reference cycle with errors.Ideally, meteorological analysis errors should be used but, as these are not available for our system during the time period of interest, a proxy was used.The proxy is based on the closest available analysis (i.e. one available 6 h prior to that used in the reference cycle), so perturbations involved a difference of analyses.The diurnal cycle was removed from the perturbations by first removing the monthly mean of each synoptic hour from the analyses.Since the resulting perturbations include the time evolution of the meteorology during 6 h as well as analysis errors, such a perturbation is expected to be larger than that due to an actual analysis error.By decomposing the departure of the CO 2 evolution in the perturbed cycle from that in the reference cycle into spherical harmonics, the spatial scales of the differences were computed.It is seen that the spatial spectra of CO 2 differences due to imperfect meteorological analyses exceeds that due to the CO 2 state itself for some scales.Thus, for spatial scales smaller than this crossover point, CO 2 is not predictable, simply due to the presence of meteorological uncertainties.The spatial scales so identified are specific to our model's domain and resolution and the monthly timescale used for averaging spectra.The spectra of the difference due to imperfect meteorology are seen to asymptote to the predictability error spectra, meaning that, for small enough spatial scales, the CO 2 field is seeing no updated information from the meteorology at all.Thus, transport error is impacted by meteorological state errors with the error being greatest for the smallest spatial scales.By comparing spectra in terms of total and zonal wavenumbers, it is seen that more information is retained in the CO 2 field in north-south direction as opposed to the zonal direction.Note that the spatial scales identified as predictable depend on the choice of fluxes used in both reference and perturbed simulations.
These experiments demonstrate that the predictability of CO 2 is limited by the presence of meteorological analysis errors as well as flux and model errors.The CO 2 observing system determines the spatial scales that can be resolved in CO 2 state estimate and in CO 2 flux estimates.However, the meteorological observing system also imposes limits on the spatial scales that can be resolved.Here these limits are found to be well below those imposed by the CO 2 observing system.However, the unresolved spatial scales are shown to increase with altitude.As the CO 2 observation density increases (for example with satellite data from GOSAT, OCO-2 and future planned missions) knowing these limits will become increasingly important.
Other components of transport error can also be isolated and compared.Considerable work has been done on identifying the contribution of model formulation errors to transport errors.Looking at the CO 2 spatial scales resolved by different model configurations provides a new tool for identifying the most important model errors.As an example, a model error was introduced by removing the transport of tracers through deep convection in a perturbed CO 2 simulation.The impact of convective tracer transport on CO 2 fields is seen to be largest in boreal spring but exceeds errors due to the use of imperfect atmospheric analyses only in the troposphere and only for wavenumbers below 5. Thus for all but wavenumbers below 5, transport errors due to imperfect meteorology are larger than those due to the lack of convective tracer transport.In other words, by comparing components of transport error, the limiting source of error for a given spatial scale can be determined.An important application of this diagnostic is as follows.With realistic flux perturbations, transport errors due to flux uncertainties could be compared to transport error due to meteorological state uncertainties.In doing so, the spatial scales of transport error dominated by flux errors and those dominated by meteorological state errors could be identified.Alternatively, if prior fluxes are used to define a reference simulation then various posterior fluxes could be used to define flux analysis increments.Specifically, the spatial scales resolved by various observing systems (e.g.GOSAT or OCO-2 vs. the surface network) could be computed and compared to the scales resolvable in the context of imperfect meteorological analyses.Indeed, such experiments are in progress and will be described in a subsequent article.
By definition, predictability experiments use a reference simulation against which perturbed simulations are compared.Thus the errors obtained here are system dependent, but if the system is representative then results are likely to be representative of other systems.Indeed, being an operational weather forecast model, GEM is routinely evaluated and we have shown that 24 h forecasts with our modified version of GEM are comparable to reanalysis products.In contrast, seasonal predictability is model dependent and is most likely related to model parameterization of subgrid-scale processes, particularly convection (Shukla and Kinter III, 2006).Subseasonal predictability is limited by the fact that most models do not capture the MJO (Waliser, 2006).Even if they did capture it, and the initial state had a strong MJO signal, forecast skill may still not be improved because of the complex interplay between MJO and other modes of variability (Lin et al., 2008).Thus, it may be useful for individual models to be able to characterize CO 2 predictability error, particularly on longer timescales, as well as the spatial scales definable in the presence of imperfect meteorological analyses.
In this work, the focus of predictability (transport) error assessment was on the random component.This is because the bias of the transport error was found to be small compared to the random component for all timescales considered here.This occurs because of the use of a control simulation when defining errors.That is, model formulation and representativeness errors are absent from such error calculations.In addition, we isolate various contributions to transport error by holding fixed other sources of error such as flux error.In reality, when comparing CO 2 model predictions to measurements all of these error sources are present and the mean error (or bias) is an important concern.Thus our methodology of isolating components of transport error provides insight into the component that is largest for a given spatial scale but may not shed light on biases that develop when all of these errors interact.
A limitation of this study is that the impact of uncertain meteorological analyses on CO 2 simulations was assessed using a 6 h shift in analysis states as a proxy for 6 h analysis errors.This was done because such estimates of meteorological analysis errors are not available but a perturbation of the approximate size and shape of a 6 h analysis error was desired.While other proxies for 6 h analysis errors could be devised, none would be any more valid.However, it is possible to directly obtain analysis and forecast errors by implementing an ensemble Kalman filter for CO 2 state estimation.Indeed, a greenhouse gas data assimilation system based on an augmented state (meteorology, constituent and fluxes) ensemble Kalman filter is now under development.This new system is called EC-CAS (ECCC Carbon Assimilation System) and also uses existing tools developed at ECCC, namely the operational global ensemble prediction system (Houtekamer et al., 2014).With EC-CAS, the atmospheric modulation of CO 2 forecast uncertainty would be directly simulated and the impact on flux estimate uncertainties could then be determined.
The Supplement related to this article is available online at doi:10.5194/acp-16-12005-2016-supplement.

Figure 1 .
Figure 1.Schematic diagram of EC-CAS forward model cycles.Meteorological analyses were precomputed by CMC's operational global deterministic prediction system using all observations collected in a 6 h window centred on the analysis time.These analyses are represented by blue boxes.CO 2 tracers are added to these analyses to form an initial condition (solid red circles) for launching a 24 h forecast (red arrows) using the coupled meteorological/tracer model starting at 00:00 UTC of each day.

Figure 2 .
Figure 2. Comparison of GEM-MACH-GHG meteorological analyses with other reanalyses for July 2009.Monthly and zonal means of differences with respect to ERA-Interim fields of GEM-MACH-GHG (left column), MERRA (middle column) and JRA55 (right column) are shown for temperature in K (top row), zonal wind in m s −1 (middle row) and meridional wind in m s −1 (bottom row).

Figure 3 .
Figure 3.Comparison of GEM-MACH-GHG with surface observations at Alert (top), East Trout Lake (middle) and Sable Island (bottom).ECCC observations (black), GEM-MACH-GHG with CT2013B fluxes (blue) and CarbonTracker-2013B (green) time series are shown for each location.

Figure 4 .
Figure 4. Mean diurnal cycle at East Trout Lake, Saskatchewan (top four panels), and Fraserdale, Ontario (bottom four panels), in 2009.Each panel shows the observed mean cycle from continuous measurements (black), CT2013B (green) and GEM-MACH-GHG with CT2013B fluxes (blue).Time is given in UTC.The four panels correspond to the months of January, July, April and October, as labelled above each panel.The grey shaded region indicates 1 standard deviation above and below observed values while the dashed lines indicate the same for the model run with the corresponding colour.

Figure 7 .
Figure 7. Time series of X CO 2 (ppm) at (a) Park Falls, USA, and (b) Wollongong, Australia, between 1 January 2009 and 1 January 2011.The dots are hourly averaged X CO 2 for TCCON (black) and smoothed GEM-MACH-GHG simulations obtained using CT2013B posterior fluxes (red).

Figure 8 .
Figure 8.Comparison of model mean profiles to NOAA aircraft observations.Observations (black curves) are from obspack_co2_1_PROTOTYPE_v1.0.4_2013-11-25 for locations over continental USA and Canada, only.Observed and modelled profiles are binned over (a) 2009, (b) March to May 2009, (c) June to August 2009 and (d) September to November 2009.CarbonTracker 2013B mole fractions (green) and GEM-MACH-GHG with CT2013B posterior fluxes (blue curves) are shown in all panels.The shaded grey regions indicate plus or minus 1 standard deviation for the observations while the dashed coloured lines indicate the same quantities but for the different model runs.Note that the x axis range differs in each panel and that ticks are every 5 ppm except for (b) where they are every 2 ppm.Sites used are Beaver Crossing, Nebraska; Bradgate, Iowa; Briggsdale, Colorado; Cape May, New Jersey; Charleston, South Carolina; Dahlen, North Dakota; East Trout Lake, Saskatchewan; Estevan Point, British Columbia; Fairchild, Wisconsin; Harvard Forest, Massachusetts; Homer, Illinois; Oglesby, Illinois; Park Falls, Wisconsin; Poker Flat, Alaska; Sinton, Texas; Southern Great Plains, Oklahoma; Trinidad Head, California; West Branch, Iowa; Worcester, Massachusetts.

Figure 9 .
Figure 9. Predictability on weather timescales during January 2009.Predictability error is defined here as the global mean zonal standard deviation of the difference in evolution of a control cycle from a climate cycle.This value is normalized by the global mean zonal standard deviation of the corresponding field in the initial condition.Predictable regimes are ones for which this ratio is much less than 1.The normalized predictability error fields for (a) CO 2 , (b) temperature, (c) vorticity and (d) divergence are plotted as a function of model vertical level (converted to approximate pressure with a reference surface pressure of 1000 hPa) and time in days since 1 January 2009.The layers labelled with white text in each panel will be used for computing the layer mean averages shown in Fig.10.

Figure 11 .
Figure11.Spectra of various fields as a function of total wavenumber.Spectra are averaged over 1 month for July 2009 (top row) and December 2009 (bottom row) and over 12 model levels.The lower and upper model levels averaged are indicated above each frame.Approximate pressure is obtained from model level by multiplying by 1000 (which corresponds to assuming a reference surface pressure of 1000 hPa).The CO 2 reference state spectra (blue curves), predictability error (black curves) and error due to a 6 h shift in analysis fields (red curves) are shown.

Figure 12 .
Figure 12.As Fig. 11 but for spectra as a function of zonal wavenumber.

Figure 13 .
Figure 13.As in 11 but for spectra in April 2009 and now errors due to the removal of convective tracer transport are shown as cyan curves.Spectra are shown as a function of total wavenumber (top row) or zonal wavenumber (bottom row).Note that only wavenumbers up to 40 are shown.
Author contributions.Saroja M. Polavarapu designed, performed and analysed the experiments and wrote the manuscript.Saroja M. Polavarapu and Michael Neish developed and diagnosed the GEM-MACH-GHG model from 2011 to 2015 and Michael Neish developed or contributed to all of the model diagnostics.Monique Tanguay implemented the tracer mass conservation scheme in GEM, while Monique Tanguay and Claude Girard devised and implemented the fixes for the continuity equation and the implementation of the dry tracer mixing ratio variable.Douglas Chan helped to analyse model results during the implementation of the dry tracer mixing ratio variable.Jean de Grandpré and Sylvie Gravel developed the GEM-MACH model on the global domain and provided input on the implementation and debugging of GEM-MACH-GHG.Kirill Semeniuk implemented the convective tracer transport scheme in GEM-MACH.Shuzhan Ren contributed to GEM-MACH-GHG model development during 2011-2014 and evaluated its PBL.Sébastien Roche and Kimberly Strong diagnosed model improvements through comparisons to TCCON measurements.

Table 1 .
TCCON stations used, their location and data reference.