Development and evaluation of CO 2 transport in MPAS-A v6.3

. Chemistry transport models (CTM) play an important role in understanding ﬂuxes and atmospheric distribution of carbon dioxide ( CO 2 ). They have been widely used for modeling CO 2 transport through forward simulations and inferring ﬂuxes through inversion systems. With the increasing availability of high resolution observations, it has been become possible to estimate CO 2 ﬂuxes at higher spatial resolution. However the computational cost of high resolution global model simulation 5 is so high that only major research and operation centers can afford it. In this paper, we implemented CO 2 transport in Model Prediction Across Scales-Atmosphere (MPAS-A). The objective is to use the variable-resolution capability of MPAS-A to enable high resolution CO 2 simulation at limited region with a global model. Treating CO 2 as an inert tracer, we implemented in MPAS-A (v6.3) the CO 2 transport processes, including advection, vertical mixing by boundary layer scheme, and convective transport. We evaluated the newly implemented model by running two sets of simulations over a 60-15 km variable-resolution 10 global domain. The ﬁrst set of simulations covers four Atmospheric Carbon and Transport-America (ACT-America) aircraft campaign seasons (2016-2018), and the simulated CO 2 is evaluated using the extensive airborne measurements from ACT. The simulation accuracy is also compared with a 27-km resolution WRF-Chem simulation and CarbonTracker (v2019) covering the same time periods. The second set of simulations covers the month of January and July of 2014, and the results are evaluated using near-surface hourly CO 2 measurements from 50 surface and tower sites across the globe. This simulation 15 accuracy is compared with European Center for Medium-Range Weather Forecasts (ECMWF) Integrated Forecasting System (IFS) global simulation conducted during the same period. Overall, the evaluation using aircraft measurements indicates that MPAS CO 2 transport model is capable of representing the observed atmospheric CO 2 structures related with the mid-latitude synoptic weather system, including the warm/cold sector distinction, boundary layer to free troposphere difference, and CO 2 enhancements along frontal boundaries. The evaluation using hourly measurements shows that the MPAS CO 2 transport model 20 is capable of achieving a same level of accuracy as the IFS 80-km resolution simulation.


Introduction
Carbon dioxide (CO 2 ) is the most important greenhouse gas, and our knowledge about its sources and sinks still have large gaps. Inversion systems are often used to infer CO 2 fluxes based on observations and chemistry transport models (CTM). Two b), Energy Exascale Earth System Model (E3SM) (Golaz et al., 2019), and Finite-Volume Cubed-Sphere model(FV3) (Putman and Lin, 2007). One benefit of local mesh refinement is enabling regional high-resolution modeling without incurring the lateral boundary condition and its associated problems, such as solution mismatches between the driving global model and the evolving regional model (Davies, 2014). 5 Model Prediction Across Scales-Atmosphere (MPAS-A) is a fully compressible non-hydrostatic global atmospheric model which uses finite-volume numeric solver discretized on centroidal Voronoi mesh with C-grid staggering of its prognostic variables (Skamarock et al., 2012;Thuburn, 2007;Ringler et al., 2010). The centroidal Voronoi mesh allows for local refinement and variable-resolution horizontal mesh which can be gradually changed from coarse to fine resolutions (Skamarock et al., 2012;Ringler et al., 2008). 10 To enable CO 2 transport modeling, we implemented in MPAS-A (v6.3) the CO 2 transport processes, including advection, vertical mixing by Planetary Boundary Layer (PBL) scheme, and convective transport. We used the newly developed model to conduct two sets of simulations over a 60-15 km variable-resolution global domain. Then the simulation results are evaluated using an extensive set of airborne observations over the eastern United States and near-surface observations from surface and 15 tower stations across the globe. The simulation accuracy of MPAS is compared with three established CO 2 modeling systems based on the same observational data: WRF-Chem (Skamarock et al., 2008;Feng et al., 2019), Carbontracker (v2019, CT2019 hereafter) (Jacobson et al., 2020) , and ECWMF IFS (Agusti-Panareda et al., 2014, 2019.
The remainder of the paper is organized as follows. Section 2 details the implementation of CO 2 transport in MPAS, Section 20 3 describes MPAS CO 2 simulation experiments and evaluation using airborne and near-surface observations, and comparison with WRF-Chem, CT2019, and IFS. Section 4 provides a summary of the model development and evaluations.
2 Implementation of CO 2 transport in MPAS This section describes the major modifications to MPAS that we made to implement CO 2 tracer transport. We represent CO 2 by its dry air mixing ratio (q co2 ) and model its atmospheric transport by adding its continuity equation in MPAS following Eq. 25 7 of Skamarock et al. (2012).
whereρ = ρ d /(∂ζ/∂z), ρ d is dry air density, ζ is the vertical coordinate, z is geometric height, t is time, and V is the velocity vector. The left hand side of the equation is the total CO 2 time tendency (∂(ρ q co2 )/∂t), and the three terms on the right hand side represent the contributions from advection, vertical mixing, and convective transport respectively. CO 2 tendency 30 from advection is modeled in flux form (Section 2.1), while tendency from vertical mixing (F bl ) and convective transport (F cu ) are modeled in uncoupled form (∂q co2 /∂t) which are coupled toρ before being added to the total tendency. Details of the three terms on the right hand side of Eq. 1 are described in the following sections. We note that because the monotonicity constraint in third-order scalar horizontal advection scheme (Skamarock and Gassmann, 2011) introduces dissipation MPAS does not use any explicit horizontal diffusion for scalar, and we did not include horizontal diffusion for CO 2 accordingly. 5

CO 2 advection
Advection is the most significant component of CO 2 atmospheric transport. Following the example of other scalars in MPAS (Skamarock and Gassmann, 2011), we model CO 2 advection as: (2) 10 where V = (u, v, w), and u, v, and w is the zonal, meridional, and vertical wind respectively. The first item on the right hand side enclosed in the square bracket is the CO 2 horizontal flux divergence, and second item is the vertical flux divergence. The horizontal flux divergence is transformed via the divergence theorem into an integral of flux over each control volume, which is modeled as: 15 where e indexes the edges of a cell and n e represents the number of edges the cell has, l e is the length of an edge, A i is the cell's areal size, F e (v H ,ρ q co2 ) is the instantaneous horizontal CO 2 flux that crosses the cell edge e, and V H = (u, v) is the horizontal wind vector. The details of MPAS instantaneous horizontal flux calculattion can found in Skamarock and Gassmann (2011). The vertical CO 2 flux divergence in Eq. 2 is calculated using finite difference 20 where F (w,ρ q co2 ) is the vertical CO 2 flux that crosses a cell's vertical face, and k indexes the vertical coordinate.

CO 2 vertical mixing
Like in WRF (Skamarock et al., 2008), a planetary boundary layer (PBL) parameterization in MPAS treats the vertical mixing of momentum and scalars not only in the boundary layer but in the entire atmospheric column. Yonsei University (YSU) PBL scheme (Hong et al., 2006) is one of the PBL schemes available in MPAS 6.3. The present YSU scheme treats vertical mixing 25 of momentum, potential temperature, and water species, but not atmospheric tracers. We modified the scheme to treat CO 2 vertical mixing.
In the YSU scheme, after the boudary layer (BL) top is determined, the vertical mixing of momentum, potential temperature, and water vapor are treated separately: above BL, local K-profile approach (Louis, 1979)) is used for vertical diffusion of momentum and scalars (Noh et al., 2003;Hong et al., 2006). Within BL, an entrainment flux at the inversion layer is included 5 for momentum and scalars diffusion. In addition, a countergradient mixing term is included for the diffusion of momentum and potential temperature to account for the convective-driven mixing (γ c of Eq. 4 in Hong et al. (2006)), but this term is not used for water vapor.
Following the treatment of water vapor, we parameterize CO 2 vertical mixing in BL as where z is the vertical distance to surface, h is BL top height, K h is vertical eddy diffusivity. Note that this formulation does not include a countergradient mixing term following the treatment of water vapor in the original YSU (Hong et al., 2006). The second term in the square bracket of Eq. 5 represents the contribution from CO 2 entrainment flux at the inversion layer, which is parameterized as: where ∆q co2 | h is the CO 2 mixing ratio difference across the inversion layer, and w e is the entrainment rate at the inversion layer calculated by Eq. A11 of Hong et al. (2006). Above BL top, vertical mixing of CO 2 is parameterized as: We use the same value for CO 2 eddy diffusivity as water vapor. The details of K h calculation can be found in the appendix 20 of Hong et al. (2006). The term ∂q co2 /∂t from Eqs. 5 is coupled with dry air density before being applied to the continuity equation (Eq. 1).

CO 2 convective transport
For convective transport, we modified the Kain-Fritsch scheme (hereafter KF) (Kain, 2004) to include the CO 2 treatment.
KF is a mass-flux convection scheme which rearranges mass in an air column using convective updrafts, downdrafts, and 25 environmental mass fluxes. Both the updraft and downdraft entrain from and detrain to the environment, thus altering the vertical profile of an air column's thermodynamic properties. We added the CO 2 convective transport as: where q co2 , q u co2 , and q d co2 are the CO 2 mixing ratio in the environment, updraft, and downdraft respectively, M u and M d are the updraft and downdraft mass respectively, ρ is the environment air density, A is the horizontal area of a cell, M = ρAδz is the mass of environmental air in a grid box, and M ud and M dd are the detrainment from the updraft and downdraft respectively.
In KF, the updraft and downdraft mass and the rates for the entrainment and detrainment are determined by a steady-state 5 plume model and a convective available potential energy (CAPE) closure assumption: 90% of the existing CAPE should be removed by the convection parameterization (Kain and Fritsch, 1990;Fritsch and Chappell, 1980;Kain, 2004). The updraft source layers are determined by a search from the model's lowest vertical level for a group of consecutive layers that is buoyant and at least 50 hPa deep (Kain, 2004). The initial value of CO 2 mixing ratio in the updraft is modeled as a pressure weighted average of the source layers: where δp k is layer's pressure depth, and q co2,k is the layer's CO 2 mixing ratio. Carbon dioxide mixing ratio of the updraft is modified by the entrainment of the environmental air through its ascent from its starting level to the cloud top.
where M ue is the updraft entrainment. The initial CO 2 mixing ratio of a downdraft (q d co2 ) is the same as that of the environment 15 (q co2 ) at the downdraft starting level and it is modified by entrainment through the downdraft descent: where M de is the downdraft entrainment.

Model evaluation
In this section we evaluate the newly developed MPAS CO 2 transport model through simulation experiments using airborne 20 and near-surface observations. After describing the simulation configuration (Sect. 3.1), we assess the model's global CO 2 mass conservation property (Sect. 3.2), then we evaluate the model accuracy by comparing MPAS simulated CO 2 with highresolution airborne measurements (Sect. 3.3) and near-surface hourly measurements (Sect. 3.4). Two sets of MPAS simulations are conducted: the first set simulations covers four ACT campaign seasons (2016-2018) for evaluation using high-resolution airborne measurements; the second set covers January and July of 2014 for evaluation using near-surface CO 2 measurements 25 and comparison with the IFS simulation results reported in Agusti-Panareda et al. (2019). In the following model evaluation, we use root mean square error (RMSE), bias (µ), and random error (STDE) as the model accuracy metrics: where o i and m i represent the observed and modeled values respectively.

Simulation experiment configuration
For all subsequent simulations, MPAS uses a 60-15km variable-resolution global mesh. Fig. 1 (Chen and Dudhia, 2001), Monin-Obukhov surface layer scheme, and WRF single-moment 6-class microphysics 15 scheme (Hong and Lim, 2006). The third-order advection is used for all scalars and CO 2 tracer. A summary of the physics parameterizations used in the simulations is given in Table 1.
Meteorology initial conditions are generated from the ERA-Interim analysis (Dee et al., 2011). To keep model meteorological fields close to the analysis, MPAS meteorological fields are re-initialized using the analysis at 00:00 UTC each day 20 throughout a simulation period. Carbon dioxide mixing ratio is kept unchanged during the meteorology re-initializations, thus a free-running simulation. This configuration is the same as that used by Agusti-Panareda et al. (2014, 2019 in their IFS global CO 2 simulations. The first CO 2 initial condition for a simulation is from CT2019 3 • × 2 • posterior dry mole fraction product and CO 2 fluxes are prepared by interpolating the CT2019 3-hourly 1 • × 1 • posterior flux product (Jacobson et al., 2020). The four fluxes from CT2019 (biosphere, ocean, fossil fuel, and fire) are interpolated to MPAS model grid and ingested at 3-hour 25 intervals throughout a simulation.

CO 2 mass conservation
For a transport model, it is very important to maintain the global CO 2 mass conservation (Agusti-Panareda et al., 2017;Polavarapu et al., 2016). To evaluate the mass conservation property of MPAS, we run a 120-hour continuous simulation and examine the change of total mass of CO 2 in the model domain. The simulation starts at 2016-08-01 00:00 UTC and ends at 2016-08-06 00:00 UTC. It uses the physics parameterizations and domain setup described in Sect. 3.1. The CO 2 dry air mixing 5 ratio field in the model is initialized using CT2019 dry mole fraction, but fluxes ingestion is turned off. Because the meteorology re-initiation will introduce dry air mass change which in turn will impact the total CO 2 mass, we run the simulation without applying the meteorology re-initialization. The impact of meteorological re-initialization and its treatment are considered in Section 3.2.2.
3.2.1 Mass conservation without meteorology re-initialization 10 Through the 120-hour simulation, model outputs are saved at 1-hour intervals using double-precision. The total mass of dry air (m air ) and CO 2 (m co2 ) is then calculated from these hourly outputs using Eqs. 15 and 16 respectively.

15
In Eqs. 15 and 16, subscript i indexes the horizontal cell, ranging from 1 to N=535, 554; subscript j indexes the vertical level, ranging from 1 to L=55. MPAS hexagon cell base areal size is A i (cell size is constant within a given column at different vertical levels), and the height of the cell is represented by h i,j (heights at different columns of a same vertical level maybe different due to the terrain-following vertical coordinate). The volume of a grid box is V i,k = A i h i,k , its dry air density (kg/m 3 ) is ρ i,j , and its CO 2 dry air mixing ratio (kg/kg) is q i,j . 20 The calculation shows that the total volume of the MPAS model domain is 1.5184682921961 × 10 19 m 3 , which is a constant through the simulation because of the model's height based vertical coordinate. At the start of the simulation, the total dry air mass is 5.053906341880670208 × 10 18 kg and total CO 2 mass is 3.079178060337270 × 10 15 kg. The total mass of dry air and CO 2 at each subsequent hour is normalized by their respective starting values and resulting variations are shown in Fig. 2. 25 The maximal variation of the total dry air mass from its starting value during the 120-hour simulation is 29, 660, 160 kg, which when divided by the total air volume represents a variation of mean dry air density approximately 1.95310 −12 kg/m 3 . The maximal variation of total CO 2 mass from its starting value is 18, 980 kg, which represents a variation in mean CO 2 mixing ratio of 1.0416 × 10 −15 kg/kg (equivalent to 6.839 × 10 −10 ppm) for a mean dry air density of 1.2 kg/m 3 .
In MPAS dry air density and CO 2 dry air mixing ratio are prognostic variables and their precision is limited by the round-off error of the model's Fortran double-precision code. Fig.2 shows that the maximal variation in dry air density is 1.95310 −12 kg/m 3 ) and maximal variation in CO 2 dry air mixing ratio is 1.0416 × 10 −15 kg/kg during the simulation period. Given that 5 the double precision float number calculation of MPAS Fortran code is capable of having 14 decimal places of precision, we consider the MPAS CO 2 transport maintains the global mass conservation to the level allowed by the machine precision.

CO 2 mass conservation during meteorology re-initialization
The last section has demonstrated that MPAS is capable of maintaining global CO 2 mass conservation. However when meteorology re-initialization is applied, dry air density fields in the model are replaced by values from the initialization files 10 generated from the ERA-Interim analysis (Dee et al., 2011). In most cases, this will cause dry air density change which in turn will introduce CO 2 mass change because CO 2 dry air mixing ratios are kept unchanged during the re-initializations.
To assess this possible change in global CO 2 mass, we conducted a 48-day MPAS simulation starting 00:00 UTC 15 July 2016 with meteorology re-initialization at 24-hour interval but without CO 2 fluxes ingestion. The global CO 2 mass is calcu-15 lated using Eq. 16 at the end of each 24-hour period which is then compared with its initial value at the start of the simulation.
The resulting global CO 2 mass variation is shown in the top panel of Fig. 3. The figure shows that during the 48-day simulation, the maximal variation of total CO 2 mass can reach approximately ±0.05%, which is 7 to 8 magnitudes larger than that without re-initialization ( Fig. 2). To restore the CO 2 mass conservation, we calculate the total CO 2 mass difference caused by a meteorology re-initialization using where notations are the same as in Eq. 16 except that ρ i,k is the dry air density after a meteorology re-initialization and ρ i,k is the value before the re-initialization. Then CO 2 mixing ratio at each grid box is scaled by where q i,k is the CO 2 mixing ratio that will be used as the initial value for the next 24-hour simulation period. 25 To test the effectiveness of this scaling method, the 48-day simulation is conducted again but with the CO 2 mixing ratio scaling after each re-initialization using Eqs. 17 and 18. The resulting variation in total CO 2 mass is plotted at the lower panel of Fig. 3. The figure shows that maximal variation is less than ±0.001%, which is small enough to be acceptable for most CO 2 transport simulations. We note that an alternative approach is to scale CO 2 mixing ratio at each grid box individually by where notation is the same as Eq. 17. This scaling approach can maintain total CO 2 mass conservation as allowed by machine precision but it will introduce artificial spatial variations in CO 2 mixing ratio. In the simulations in the following sections, we chose to use the first scaling approach to avoid the artificial CO 2 mixing ratio variation by accepting the small change in total 5 CO 2 mass.

Model valuation using airborne measurements
In this section, we evaluate the MPAS CO 2 simulation accuracy using an extensive high resolution CO 2 observation data acquired through the Atmospheric Carbon and Transport-America project (ACT-America, henceforth, referred to as ACT). ACT is a National Aeronautics and Space Administration (NASA) Earth Venture Suborbital 2 (EVS-2) mission, and its goal is to 10 improve atmospheric inversion estimates of CO 2 and CH 4 through extensive airborne measurements over the eastern United

Overall accuracy of MPAS modeled CO 2
We use the ACT airborne measurements to evaluate MPAS CO 2 simulation regarding its overall accuracy and its performance measured by three model evaluation metrics proposed by Pal et al. (2020). To provide an objective reference, we also compare 5 MPAS performance with two established CO 2 model systems: WRF-Chem (Skamarock et al., 2008) and CT2019 (Jacobson et al., 2020) using the same set of airborne measurements.
WRF-Chem is an online chemistry transport model based on the regional NWP WRF ( 20 We use the ACT 5-second averaged CO 2 measurement dataset (Davis et al., 2018b) , which has a horizontal resolution approximately 500 m given the average aircraft velocity. MPAS simulated CO 2 fields are interpolated in time and space to match each 5-second airborne data points. WRF-Chem simulated CO 2 fields are also interpolated to match the ACT 5-second data point using the same approach as MPAS. CarbonTrack CO 2 used for the evaluation is obtained from CarbonTrack ObsPack ((Masarie et al., 2014)), which is the CT2019 posterior mole fraction interpolated to the ACT 5-second data points. one for the boundary layer (BL) and another for he free troposphere (FT). RMSE and bias of the modeled CO 2 from the three models (MPAS, WRF-Chem, and CT2019) are then calculated for each flight day, with the calculations for the BL and FT carried out separately. In addition to the RMSE and bias, standard deviation is also calculated for the observed and modeled CO 2 for assessing each model's representation of CO 2 spatial variability.
Using RMSE as an accuracy metric, we compare MPAS CO 2 simulation with WRF-Chem and CT2019 in Fig. 4. The figure shows that all three models have higher accuracy (lower RMSE) in FT than BL, which is most likely because flux errors impact simulated CO 2 in BL more than in FT. The figure also shows that compared with WRF-Chem, MPAS has a similar magni-5 tude of RMSE overall in both BL and FT, with some exception where substantial differences exist between the two models.
Averaged over the 97 flight days, the mean RMSE is 4.49/1.85 ppm (BL/FT) for MPAS, and 4.91/1.71 ppm for WRF-Chem, indicating MPAS achieved a slightly higher accuracy than WRF-Chem in BL which maybe partially because the former has a higher horizontal resolution (15 km) than the latter (27 km). In free troposphere, WRF-Chem resulted in a slightly higher accuracy than MPAS maybe because that it applied meteorological nudging while MPAS did not. Compared with CT2019, MPAS 10 simulations resulted in larger RMSE in the majority of flight days, and the differences in RMSE between the two models are large in BL than FT. Averaged over the 97 flight days, the mean RMSE for CT2019 is 3.36 ppm and 1.42 ppm in BL and FT respectively, which are substantially lower than MPAS (4.49 ppm/1.85 ppm in BL/FT).
Next we examine how well each model represents the CO 2 spatial variability (as measured by standard deviation). The above model evaluation using airborne observations shows that the MPAS CO 2 transport model is capable of achieving a similar level of accuracy as 27 km resolution WRF-Chem simulation. Compared with the global inversion system CT2019, although MPAS resulted in higher RMSEs, it achieved better estimations of the observed CO 2 spatial variability.

Model representation of CO 2 difference between warm and cold sectors
Through analyzing the ACT Summer 2016 campaign data, Pal et al. (2020) identified three consistent features in CO 2 mole fraction and proposed to use these features as transport model assessment metrics. The three features are the differences between the warm and cold sectors, the difference between the boundary layer and free troposphere, and the CO 2 enhancement 30 bands in the vicinity of frontal boundaries. Here and in next two sections, we evaluate how MPAS simulated CO 2 represents the three features.
Using the ACT maneuver flag dataset (Pal et al., 2020), we identified flights that crossed a weather front and their associated warm and cold sectors. The CO 2 mole fraction statistics for the warm and cold sectors are calculated from the aircraft measurements and the modeled CO 2 by MPAS, WRF-Chem, and CT2019, respectively. The results are shown in Fig. 6, which summarizes the statistics of CO 2 mole fraction differences between the warm and cold sectors measured by 15 front-crossing flights: 10 from the summer 2016 season and 5 from the winter 2017 season. The figure confirms that warm sector has higher 5 average CO 2 mole fraction in the boundary layer than the cold sector during summer 2016 as reported by Pal et al. (2020).
The figure also shows that the average CO 2 mole fraction in the warm sectors are lower than than the colder sectors in winter 2017, opposite to the summer 2016.  model is capable of well representing the observed CO 2 difference between the warm and cold sectors, and its accuracy in this respect is similar to WRF-Chem and CT2019.

Model representation of CO 2 vertical difference
The second feature identified by Pal et al. (2020) is the vertical difference of CO 2 mole fraction between BL and FT. Dur-20 ing ACT campaign season, two research aircraft (B200 and C130) took many vertical profile measurements during take off, landing, spiral up and down, and inline ascend and descend maneuvers (Pal, 2019). These profile observations characterize the vertical variation of atmospheric CO 2 mole fraction. From the vertical profile measurements taken during the summer 2016 season, Pal et al. (2020) calculated the mean CO 2 mole fraction in the boundary layer (BL) and free toposphere (FT), denoted

25
They found that ∆[CO 2 ] tend to be positive in the warm sector and negative in the cold sector. In this section, we evaluate how well MPAS represent the BL-to-FT CO 2 difference and compare its performance with WRF-Chem and CT2019.
Using the ACT maneuver flag dataset (Pal et al., 2020), we identified all vertical profiles taken during the four campaign seasons, from which we selected profiles that meet two criteria: (1) a vertical profile must include at least 20 5-second mea-  Table 5 which shows that: MPAS has lower RMSEs than WRF-Chem in all but the summer 2016 season, and 20 it has lower RMSEs than CT2019 for the winter and fall of 2017 but higher in the other two seasons. Regarding the standard deviations, MPAS compares with the observation better than WRF-Chem in all four season, and better than CT2019 in all but the summer 2016 season. In summary, the above model evaluation and comparison demonstrate that MPAS CO 2 transport model is capable of representing the aircraft observed CO 2 difference between boundary layer and free troposphere at least as accurately as WRF-Chem and CT2019. To provide an quantitative comparison between the three models' accuracy regarding the frontal boundary CO 2 variation, we calculated RMSE and standard deviation for each of the 48 level-leg flights shown in Fig. 9. For each flight, RMSE is calcu-25 lated for each of the three models as compared with the aircraft observations, and standard deviation is calculated for both the models and aircraft observations. The resulting RMSEs from MPAS are compared with WRF-Chem and CT2019 in Fig. 10.
The figure indicates that MPAS RMSEs in general are similar in magnitude to CT2019 and lower than WRF-Chem. Averaged over the 48 cases, the mean RMSE is 4.63 ppm for MPAS, 4.68 ppm for CT2019, and 5.86 ppm for WRF-Chem, indicating MPAS perform as well as the other two models as measured by the RMSEs. We also assess how well the three models estimate 30 the spatial variability of the frontal boundary CO 2 by comparing the standard deviation (σ) of the three models with the observations (Fig. 11). The figure shows that MPAS and WRF-Chem represent the spatial variability in the 48 level-leg flights better than CT2019 which substantially underestimates the variability in the majority of the cases. As horizontal resolution impacts a model's ability to represent small scale spatial variability (Agusti-Panareda et al., 2019), the coarser resolution of CT2019 (1 • × 1 • over the North America) is likely the primary cause of its underestimation of the frontal boundary CO 2 variability.
3.4 Model evaluation using near-surface hourly CO 2 observations The MPAS CO 2 transport model has been evaluated using the extensive high-resolution aircraft measurements from four ACT campaign seasons. Considering that these aircraft measurements were all acquired during the day time hours over the eastern   Table S1 of the supplement material. The table shows that RMSEs of the MPAS simulated hourly CO 2 ranges from 0.36 ppm at the spo to 25.96 ppm at ssl station. In comparison, the IFS simulations also resulted in a much lower RMSE at the spo than ssl, the latter of which has a RMSE of 5.83 ppm from 9 km resolution sim-  Table S2 (supplement 10 material) and Fig. 13. Comparing with the January simulation, Fig. 13 shows that both MPAS and IFS have larger RMSEs in July, suggesting the impacts of large uncertainty in biospheric carbon fluxes during boreal summer. Fig. 13 also shows that the MPAS simulation has larger RMSE and STDE than the 9 km IFS simulation but similar to the 80 km IFS simulation, a similar pattern as in January. re-initialization and fluxes, and the results show that MPAS is capable of maintaining CO 2 mass conservation to the limit of machine precision. A 48-days simulation with meteorology re-initialization indicates that changes in dry air density during the re-initialization causes changes in global total CO 2 mass, and a scaling method applied after each re-initialization is able to reduce the global CO 2 mass variation to less than ±0.001% of its initial value.

30
The accuracy of MPAS CO 2 transport is evaluated using the extensive high-resolution aircraft measurements from four ACT campaign seasons during 2016-2018. Compared with a 27 km resolution WRF-Chem simulation and CT2019 posterior CO 2 mole fraction, MPAS simulated CO 2 achieves a comparable level of accuracy (as measured by RMSE). Further evaluation using three metrics proposed by Pal et al. (2020) shows that MPAS simulation is capable of representing the observed CO 2 features as accurately as the WRF-Chem simulation and CT2019.
A second set of MPAS simulations for the month of January and July of 2014 are evaluated using near-surface hourly CO 2 measurements from surface and tower stations across the global. The resulting statistic, including RMSE, random error, and The model evaluations using the airborne and near-surface measurements, indicates that the newly developed MPAS CO 2 10 transport model is capable of achieving a comparable level of accuracy with the more established CO 2 modeling systems, including the regional model system WRF-Chem, the operational assimilation system CarbonTracker, and the lower resolution (80 km) simulation of ECWMF IFS global CO 2 modeling sytem. Although further improvements are expected, the MPAS CO 2 transport model has the potential to contribute to improving our knowledge of the atmospheric CO 2 transport and fluxes.

Code and data availability 15
Source code for MPAS-A CO 2 transport model v6.3 can be retrieved at https://doi.org/10.5281/zenodo.3976320 (Zheng, 2020      27 https://doi.org/10.5194/gmd-2020-265 Preprint. Discussion started: 7 October 2020 c Author(s) 2020. CC BY 4.0 License. qq q q q qq q q q q qq q q q q qqq q qq q q q q q qq qqq q q q qqqq q qq q q qq q qqqq q qqqq q q q qqq qq qq q qq q q q q q q q q q qq q q q q q qqq q qqqq q q q q qq q q q q q q q q q qq q q q qqq q 0 20 40 60 80 100 120 0e+00 4e−12 8e−12 (a) Total dry air mass change as ratio to the initial total dry air mass qq q q q q q q qqqqq q q q q qq q qqq qq q qq q q qq qq q q q q q q q q qqq q q qqq q q q qq q q q q qq qq qqqqqq q q qqq q q q q qq qqq q q q q q qq q q q q q q q qq qq q q q q qq q q qq qq q qq q q 0 20 40 60 80 100 120 0e+00 4e−12 8e−12 Index (co2 − co2_init)/co2_init (b) Total CO2 mass change as ratio to the initial total CO2 air mass Number of hours since the start of simulation Figure 2. Variation of total dry air mass (a) and total CO2 mass (b) as the ratio to their respective starting values during a 120-hour continuous MPAS simulation without meteorology re-initialization. The X-axis represents the number of hours after the start of the simulation, and the Y-axis the ratio of the total mass change to the starting values. (a) Without CO2 mixing ratio scaling q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q  . Variation of total CO2 mass as the ratio to its starting value during a 48-day MPAS simulation with meteorology re-initializations at 24-hour intervals. The top figure is from the simulation without applying CO2 mixing ratio scaling as described in Sect. 3.2.2, and the bottom figure is from the simulation with the scaling. In each figure, X-axis represents the number of days after the start of the simulation, and Y-axis the ratio of total CO2 mass variation to its starting value.