Mass-conserving coupling of total column CO 2 ( X CO 2 ) from global to mesoscale models : Case study with CMS-Flux inversion system and WRF-Chem ( v 3 . 6 . 1 )

Quantifying the uncertainty of inversion-derived fluxes and attributing the uncertainty to errors in either flux or transport continue to be challenges in the characterization of surface sources and sinks of carbon dioxide (CO2). It is also not clear if fluxes inferred in a coarse-resolution global system will remain optimal in a higher-resolution modeling environment. Here we present an off-line coupling of the mesoscale Weather Research and Forecasting (WRF) model to optimized biogenic CO2 fluxes and mole fractions from the global Carbon Monitoring System inversion system (CMS-Flux). The coupling framework 5 consists of methods to constrain the mass of CO2 introduced into WRF, effectively nesting our North American domain within the global model. We test the coupling by simulating Greenhouse gases Observing SATellite (GOSAT) column-averaged dryair mole fractions (XCO2) over North American for 2010. We find mean model-model differences in summer of ⇠0.12 ppm. While 85% of the XCO2 values are due to long-range transport from outside our North American domain, most of the modelmodel differences appear to be due to transport differences in the fraction of the troposphere below 850 hPa. The framework 10 methods can be used to couple other global model inversion results to WRF for further study using different boundary layer and transport parameterizations.

observations.The CMS-Flux inversion system uses the forward GEOS-Chem global chemical transport model (Nasser et al., 2010(Nasser et al., , 2011) ) driven by meteorological fields from the NASA Goddard Earth Observing System, Version 5 (GEOS-5) data assimilation system (Rienecker et al., 2008).CO 2 is simulated as a passive tracer forced by emissions from fossil fuel, biomass burning, shipping and aviation, biogenic land, and ocean surface fluxes.Nasser et al. (2011) used GEOS-Chem as the transport model in a Bayesian synthesis inversion for surface CO 2 fluxes constrained by mid-tropospheric CO 2 from the Tropospheric Emission Spectrometer.The GEOS-Chem adjoint model (Henze et al., 2007) optimizes the biogenic surface fluxes and atmospheric CO 2 mole fractions to be consistent with XCO 2 from a satellite observing system (Liu et al., 2014(Liu et al., , 2015)).The GEOS-Chem adjoint has been used to estimate carbon monoxide emissions (Kopacz et al., 2009(Kopacz et al., , 2010) ) and to attribute direct ozone radiative forcing (Bowman and Henze, 2012).Here we use the imposed surface fluxes and the optimized biogenic fluxes and atmospheric CO 2 mole fractions from a coarse resolution global CMS-Flux inversion (horizontal resolution: 4 latitude x 5 longitude; vertical resolution: 47 levels to 0.01 hPa) that assimilated GOSAT XCO 2 for 2010.Details of the imposed and prior fluxes can be found in Liu et al. (2014) and are summarized for North America in Table 1. Figure S1 shows a schematic of the CMS-Flux system and the WRF-CMS interface for our experimental setup.

Weather Research and Forecasting with Chemistry Model (WRF-Chem)
The mesoscale model is WRF-Chem v3.6.1 (Powers et al., 2017;Skamarock et al., 2005;Grell et al., 2005;Fast et al., 2006) with the modification described in Lauvaux et al. (2012) to transport greenhouse gases as passive tracers.Trace gas boundary conditions are provided from a global model at 6-hourly intervals and surface fluxes introduced hourly.This work will describe (Section 2.3 below) the framework for introducing the boundary conditions from the global model into the WRF domain in a manner designed to preserve the vertical distribution of CO 2 at the WRF domain edges.Lauvaux et al. (2012) used WRF-Chem for the forward transport model in an inversion for CO 2 fluxes in the USA Upper Midwest region as part of the Mid Continent Intensive (MCI) campaign.This WRF-Chem implementation was also used in Lauvaux and Davis (2014), who investigated the impact of introducing column-averaged XCO 2 into regional inversions that typically use only CO 2 observations measured in the planetary boundary layer.
In the experiment described here, the regional domain contains most of North America in a Lambert Conformal projection at 27 km horizontal resolution.The model has 50 levels up to 50 hPa with 20 levels in the lowest 1 km.The model meteorology is initialized every 5 days with 6-hourly ERA-Interim (Dee et al., 2011) reanalysis at T255 horizontal resolution (⇠80 km resolution over North America).Each meteorological re-initialization is started at a 12-hour setback from the end of the previous 5-day run.The first twelve hours of the new run are then discarded.We also employ the alternative lake surface temperature initialization as described in the WRF User Guide (Weather Research & Forecasting ARW Version 3 Modeling System User's Guide, January 2015, pp3-26).Choices of the model physics parameterizations used in this experiment are documented in Table 2.
Carbon dioxide from surface fluxes from the CMS-Flux inversion (optimized biogenic fluxes and imposed surface fluxes from the ocean and emissions of fossil fuel, biomass burning, and ship bunker fuel; Table 1) are carried as individual tracers in WRF-Chem.Background CO 2 mole fractions, supplied as boundary conditions from the CMS-Flux optimized mole fractions, Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-342Manuscript under review for journal Geosci.Model Dev. Discussion started: 1 February 2019 c Author(s) 2019.CC BY 4.0 License.are in a separate tracer.For analyses requiring total CO 2 or XCO 2 , the surface flux tracers are summed, dried of water vapor content, and added to the boundary condition tracer.This multiple-tracer strategy allows for inspection of the separate tracers throughout the model integration.If total column values are required, the region above the 50 hPa top of the WRF model is populated with the appropriate value from 50 hPa to 0.01 hPa from the global model.In this experiment, we do not include the non-surface fluxes of CO 2 from the CMS-Flux inversion (aircraft source and chemical source; Table 1), but these are included in the boundary condition mole fractions from the global domain.A test, reducing these non-surface sources to surface emissions in the WRF domain, made, at most, a 0.03 ppm difference in the simulations of GOSAT XCO 2 .
The final preparation step for our experiment is the population of the WRF atmosphere with CO 2 from the global model.
It takes approximately a month's integration in model time to distribute the contributions from the boundaries throughout the WRF atmosphere.For this experiment, we started the model integration in early December 2009, but begin all analyses in January 2010.

Mass Conserving Coupling Framework
In this experiment, we wish to compare simulated column-averaged CO 2 from both the WRF execution and the CMS-Flux products with the ACOS GOSATv3.5 soundings (O'Dell et al., 2012) for North America.This requires that we ensure that the mass of CO 2 introduced into the WRF domain from the CMS-Flux GEOS-Chem global model be conserved, both at the surface (fluxes) and at the boundaries (atmospheric mole fractions) of the WRF model domain.The challenges for mass conservation include differences in model horizontal grid resolution, implied grid surface elevation, and vertical grid discretization.The strategy described here can also be used to couple the WRF regional model to other global models with a minimal amount of customization specific to the global model.

Mass Conservation of Surface Fluxes
We apply domain-wide scaling factors for each surface flux to conserve the mass of CO 2 introduced into the domain by the surface fluxes.First, we create a map projection for translation from the CMS-Flux GEOS-Chem 4 x 5 grid to the WRF 27 km grid (Figure S2 illustrates the domain extent of the WRF Lambert Conformal projection).After assignment of WRF grid cells to the global grid, we compute monthly scaling factors for each surface flux as follows: 1. Calculate the sum of the mass exchange in the global model grid cells assigned to the WRF domain.
2. Calculate an initial domain-wide sum of mass exchange for the WRF grid cells using the assigned global model grid cells.
3. Compute a domain-wide scaling factor as the ratio of the results of step 1 and step 2.

Mass Conservation at the Domain Boundaries
The challenge in the case of the domain boundaries is not only the difference in horizontal resolution and grid type, but also the difference in vertical discretization schemes.WRF uses a terrain-following, hydrostatic-pressure vertical eta coordinate with a fixed model top.The CO 2 boundary mole fractions for our experiment are from a GEOS-Chem model with 47 hybrid-sigma layers corresponding to the GEOS-5 (MERRA) reduced vertical grid (Rienecker et al., 2008).This grid is terrain-following from the surface up to ⇠170 hPa, with fixed pressure levels from 170 hPa to the top of the model.Depending on the surface pressure and terrain, there may be 9 layers below 1 km, compared to 20 layers below 1 km for our WRF domain.To approximate the vertical distribution of CO 2 in the global model source in the WRF grid, we follow these steps for each WRF boundary grid cell: 1. Compute an equivalent global model surface pressure for each WRF boundary grid cell using standard bi-linear interpolation (http://numerical.recipes;Interpolating in two or more dimensions) using the four global model surface grid cells whose centers are closest in latitude and longitude to the WRF grid cell center.
2. Compute the equivalent global model pressure column using this derived surface pressure and the standard vertical grid discretization for the global model, in this case the GEOS-5 hybrid sigma-pressure ap and bp parameters and algorithm for the 47-layer reduced vertical grid (http://wiki.seas.harvard.edu/geos-chem/index.php/GEOS-Chem_vertical_grids).
3. Compute an equivalent WRF pressure column using this derived surface pressure and the znu WRF vertical resolution discretization vector.
4. Assign a source global model level for each receiving WRF model level in the WRF boundary grid (Figure 1).This is determined by the global model level edge pressures between which the derived WRF midpoint layer pressure falls.This computation is done in log space with the respective pressure columns in units of Pa. 5. Use simple linear interpolation (if necessary) between global model levels to smooth out a poorly mixed flux signal.We used this technique for CMS-Flux GEOS-Chem where the first four or five model levels in WRF are sourced from the first GEOS-Chem level; this source layer often shows the immediate result of the surface flux as distinct from several well-mixed layers above the surface layer.6.The result of this transfer of vertical CO 2 columns from the global model to the WRF boundary grid cells is introduced into WRF via the WRFCHEMBC functionality used for meteorological boundary conditions, similar to the approach of (Ahmadov et al., 2007(Ahmadov et al., , 2009)).
We use the GEOS-Chem surface pressure for two reasons.First, it is the mass in the GEOS-Chem column that we want to introduce into the WRF model domain, and, second, it is not possible to match the surface pressures between the two models due to different horizontal grid resolutions, model grid surface elevations, and driver meteorology.To test the adequacy of this method, we compute pressure-weighted column-averaged XCO 2 along the WRF boundaries, independently compute the same quantity from the global model up to 50 hPa, and compare the results.An example of this comparison for a day in early June 2010 is shown in Figure 2. The surface layer of the western and eastern boundaries of the WRF domain is predominantly ocean, where we do not expect significant model grid elevation or surface pressure differences.This is evident for this example date in Figure 2b.On the other hand, we do expect some differences in the southern and northern boundaries, particularly in mountainous areas (Figure 2a).Our algorithm produces model-model differences of column-averaged CO 2 at the boundaries of less than 0.1 ppm on most of the boundaries and less than 0.3 ppm in the high-terrain regions of the northern boundary.

Model Comparison to Observations
While a primary goal of this model coupling experiment is to compare WRF and CMS-Flux model-simulated XCO 2 GOSAT soundings, it is comparisons to other observations that provide verification and insight into model behavior.In addition to the GOSAT XCO 2 , we compare to XCO 2 from a Total Carbon Column Observing Network (TCCON) site for times of GOSAT overpasses.We also compare model meteorological winds to observations from selected North American rawinsonde sites.

GOSAT XCO 2
We simulate in both the CMS-Flux and WRF model atmospheres more than 13,000 good quality ACOSv3.5 LITE GOSAT soundings (http://co2.jpl.nasa.gov)within the WRF North American domain during 2010.The distribution of these soundings varies in space and time, with no coverage in ocean areas or north of 60 N in winter.We follow the same sampling scheme in both WRF and CMS-Flux model CO 2 mole fractions: 1. Locate the nearest model grid cell in time and space to the GOSAT sounding.
2. Isolate the column of CO 2 at that grid cell, calculate the dry air mole fractions.
3. Interpolate the simulated CO 2 to the 20 pressure levels specified in the sounding's ACOSv3.5 averaging kernel algorithm and a priori profile (O'Dell et al., 2012;Rodgers and Connor, 2003) as shown in the ACOS 3.5 User Guide.
4. Apply the averaging kernel to the interpolated profile for the full column XCO 2 or use the a priori profile pressure-level weights for partial column analysis.
Using the pressure-level weights versus the complete averaging kernel for the whole column yields results within a few tenths of a ppm for the full columns in our model atmospheres.For a priori profile pressure levels above 50 hPa, we use the CMS-Flux

Lamont TCCON XCO 2
The Total Carbon Column Observing Network (TCCON; Wunch et al. ( 2010)) provides measurements of XCO 2 from the earth's surface at a global network of sites.Two TCCON sites within our North American domain were operating in 2010.
The Park Falls site has significant drop-outs, especially in summer, so we choose the Lamont TCCON site (36.6 N, 97.49W) for comparison (Wennberg et al., 2017).We average GOSAT soundings in a box of 6 latitude and 12 longitude centered at the Lamont site and match them to the TCCON data averaged for the hour of the GOSAT overpass of this regional box.We also computed hourly-averaged simulated XCO 2 from the models for these GOSAT soundings to compare to both GOSAT and TCCON XCO 2 .During our period of comparison, summer of 2010, this region covers a large gradient of surface CO 2 fluxes.Our goal is not to evaluate the GOSAT XCO 2 relative to TCCON observations, but rather to illustrate how well the model-simulated XCO 2 follow the general tendencies of the observations.More rigorous evaluation of the GOSAT XCO 2 would require use of coincidence criteria as shown in Wunch et al. (2011) or Basu et al. (2013).

Horizontal Wind at Selected Rawinsonde Sites
We also make use of the NOAA archive of rawinsonde data (http://www.esrl.noaa.gov/raobs/fsl-format-new.cgi; Schwartz and Govett (1992)) for model comparisons to observed wind speed and direction at mandatory reporting levels.The goal in this case is to identify possible causes of transport differences between models and observations.Comparisons of 00 UTC soundings (16)(17)(18)(19) LST, depending on the site) are summarized to show annual and seasonal biases.

Evaluation of Transport Differences in Model-Simulated GOSAT XCO 2 Soundings
We compare the CMS-Flux and WRF XCO 2 model simulations with each other in Figure 3 and with the GOSAT soundings in Table 3.If we have reached our target of conserving CO 2 mass entering the WRF domain from the CMS-Flux optimized fluxes and boundary conditions, then we expect the model-simulated XCO 2 to be similar.Model-model differences can be related to transport (horizontal transit times from different driver meteorology and vertical mixing from boundary layer processes) or model resolution (heterogeneous surface characteristics).The seasonal mean model-model differences in XCO 2 simulations are largest in winter (mean, -0.30 ppm; RMSD, 0.51 ppm) and smallest in summer (mean, 0.12 ppm; RMSD, 0.72 ppm).
These seasonal distributions of differences are all sharply peaked around zero with a few large outliers (>5 ppm) in summer.
Comparisons of model simulations to the GOSAT soundings are summarized in Table 3.In spite of a few extreme outliers in winter and spring, the inner quartile ranges (IQRs) for both models relative to GOSAT are approximately 2 ppm, with nearly the same seasonal mean differences (CMS-Flux range [0.00, 0.31] ppm; WRF range [-0.04, 0.38] ppm).Median differences are GOSAT soundings assigned to a single grid cell in GEOS-Chem are represented by many grid cells in the higher resolution WRF grid.However, we do not see any differences at this summary level.There are six individual GOSAT XCO 2 soundings with model-GOSAT differences greater than 10 ppm; these are all in challenging terrain in the western United States.
These general comparisons, however, do not show differences in spatial coverage by season of the GOSAT XCO 2 soundings or any spatio-temporal variations in model-model residuals.To address this, we aggregate the model-model differences to the CMS-Flux GEOS-Chem 4 x 5 grid to map the mean seasonal differences in each grid cell (Figure 4).We report in Figure 4 only those grid cells with more than 10 GOSAT soundings during each season shown.There are consistent small negative residuals in the WRF XCO 2 simulations relative to CMS-Flux in most of the continent in winter and in the south and east in fall.The WRF -CMS-Flux differences are slightly positive in the northwest in fall.The pattern in spring is mixed.There are strong positive differences in the Pacific Northwest in summer.The underlying CMS-Flux optimized biogenic flux shows a very strong source of CO 2 in the upper Pacific Northwest, and a correspondingly strong sink in the MidWest and East in July.
These sources and sinks are evident in the WRF model layers closest to the surface, but not always in the CMS-Flux optimized mole fractions in the surface layers in the same locations.On some days, the surface layer flux contributions in the CMS-Flux mole fractions are in grid cells adjacent to the emitting grid cells, suggesting some model-model differences in mixing and transport within the boundary layer.There are also sharp gradients in terrain height in this region that may contribute to the model differences, although this does not appear to be an issue in the Rocky Mountain region of the US West.In general, there is more variability in the spatial differences during the growing season.These maps also show how the spatial distribution of the GOSAT soundings changes by season.Note that there are no ocean soundings and no soundings in the high latitudes in winter.Examination of variances of the model-simulated soundings on this spatial scale show no clear differences between models, other than that there is less spatial variance in the models than in the GOSAT soundings.
Having documented the model-model differences in column-averaged XCO 2 at seasonal scales, we next compare the distribution within the columns of the CMS-Flux and WRF XCO 2 simulations.We expect to find the most differences in the active growing season, so focus the comparison on summer.We divide the columns into upper and lower portions, with the lower column corresponding to the three pressure levels closest to the surface in the ACOS algorithm, and with the upper column consisting of the remaining 17 levels.The ACOS averaging kernel pressure level weights specify ⇠12.5% of the columnaveraged value from the three pressure levels closest to the surface.The elevation above the surface of this split varies with terrain, but roughly corresponds to 850 hPa.As justification for this division, we highlight an example CO 2 profile at the location of the LEF tall tower in Park Falls, Wisconsin, USA (45.95 N, 90.27 W) at the time of a GOSAT overpass on 27 August, 2010 (Figure 5).In this example, the XCO 2 simulations from both the WRF total CO 2 and the CMS-Flux optimized mole fractions agree with each other and with the ACOS GOSAT XCO 2 (Figure 5a).The GOSAT XCO 2 value is 385.526ppm with an uncertainty of 1.075 ppm; WRF and CMS-Flux simulated values are 385.706and 385.245 ppm, respectively.We decompose the WRF total column CO 2 into contributions from the global model (light blue boundary conditions profile in Figure 5a) and the contributions of the flux tracers (Figure 5b).This is a location far from the boundaries of the WRF regional domain, so the boundary conditions profile is the result of long-range transport, not recent inflow.At this location and time, the flux contributions to the total CO 2 profile are predominantly below 850 hPa.The minor fluxes (ocean, and combined biomass burning, biofuel, and ship bunker fuel emissions) contribute very little to the overall column value.The biogenic flux and the imposed fossil fuel emissions constitute almost all the flux portion of the column value.The transported boundary conditions account for more than 85% of the column-averaged CO 2 .For applications using total column CO 2 , it is important that this coupling of regional to global model be done correctly.
We illustrate the model-model differences for the simulated GOSAT XCO 2 in summer using this ⇠850 hPa division in Figure 6.Differences shown are for sub-column-averaged CO 2 within the lower and upper portions of the column, computed as the sum of the interpolated mole fractions at each ACOS pressure level times the ACOS-provided weight at that pressure level.The lower portion of the column accounts for ⇠50 ppm of the total column-averaged CO 2 .Comparing spatial patterns of differences of the sub-columns versus the total column in Figure 4c, we see that the WRF positive residual in the Pacific Northwest is the result of larger mole fractions in both upper and lower sub-columns.This is an area with a relatively low count of GOSAT soundings (Figure S3) and very large biogenic source flux.In the eastern half of North America, there is a dipole effect, with WRF simulations having lower mole fractions in the lower sub-column and higher values in the upper sub-column compared to the CMS-Flux simulations.There are non-homogeneous patches of strong uptake in the biogenic flux in eastern North America in July, generally matching the spatial pattern in the lower sub-column.The upper sub-column pattern more closely resembles the total column pattern in Figure 4c.Model-model differences in other seasons are minimal in both parts of the column.The summer patterns illustrate that agreement of full column XCO 2 simulations does not necessarily imply that the distribution of the CO 2 is the same in the two models.This further suggests model-model differences in transport within the columns, possibly due to more vertical mixing in the CMS-Flux GEOS-Chem model.The CMS-Flux vertical profile in Figure 5a shows an example of this behavior.

Temporal Evaluation of Model-Simulated XCO 2 at the Lamont TCCON Site
An independent comparison of column-averaged XCO 2 for GOSAT and the model simulations is possible using the XCO 2 observations at the Lamont TCCON site.GOSAT passes near the Lamont TCCON site around 19 UTC.We select the GOSAT soundings near Lamont, average the XCO 2 values by day and report them along with the 19 UTC hourly average of the TCCON observations.We present the weighted averages and standard deviations of both sets of observations in Figure 7a.The error bars in Figure 7a represent the root mean square uncertainties of the selected TCCON and GOSAT soundings included in the hourly averages.In the example shown for days in July and August 2010 when both observations are available, there are 1-21 good GOSAT soundings and 3-34 TCCON XCO 2 observations during the 19 UTC overpass.Despite the presence of outliers in the early summer , daily average GOSAT and TCCON XCO 2 soundings agree within 1-2 ppm, with discrepancies likely due to sampling and representation errors.The mean residuals or ACOS GOSAT, CMS-Flux, and WRF XCO 2 relative to TCCON for this period are 0.088 ± 1.856 ppm, 0.868 ± 1.042 ppm, and 0.965 ± 0.981 ppm, respectively (Figure 7b).The residuals from the model simulations are generally, but not exclusively, positive and are more similar to each other than to either of the observations.This implies that model transport differences are small compared to other model-data differences at this time and location.The synoptic-scale variability in the summer atmospheric CO 2 mole fractions is well represented during short-term events (e.g., near DOY 218).Despite the coarser resolution of the CMS-Flux mole fractions compared to WRF, both models are able to capture the trend and variability observed by the Lamont TCCON site.

Horizontal Mean Wind
We selected a group of rawinsonde sites (Table 4; Figure 8) including west coast, east coast, Gulf Coast, and mid-continent North American sites for comparison of model and observed wind speed and wind direction at mandatory rawinsonde reporting levels.Our intent is to identify any regional distinctions in biases or variability that might help to explain the CO 2 model-data differences.Each of these sites had more than 11 months of data for 00 UTC soundings in 2010 at the mandatory reporting levels of 925 hPa, 850 hPa, 500 hPa and 250 hPa which we used in our analysis.Continental mountain and high plains sites were not included because they lack 925 hPa data, and cannot easily be compared with mandatory levels at lower elevation sites.
Model locations were assigned as the nearest grid cell with similar surface height.For some coastal sites, the representative model sites were moved one grid cell toward the ocean to achieve this.In Figure 9, we show the annual mean wind speed bias at these sites from the WRF forward run (initialized with the ERA-Interim reanalysis) and the GEOS-5 reanalysis used by the CMS-Flux system.Recall that although the WRF run is initialized with reanalysis data, it is free-running during each 5-day run, and so may not conform to the reanalysis.Wind speed bias, RMSE, variance ratio and correlation skill score ( Cov(model,obs)  p V ar(model)V ar(obs) ); von Storch and Zwiers (1999)) are summarized in Table 5 for mandatory reporting level 925 hPa.(See Tables S1-S3 for results at other reporting levels.)In general, coastal sites are difficult to simulate in both models, with seasonal differences that tend to cancel each other in WRF but not in the GEOS-5 meteorology.With the exception of the coastal sites, WRF wind speed error is positively biased at the lower levels, but this bias largely disappears with height.Most notable is that GEOS-5 underestimates the wind speed for nearly all of this selection of sites at all levels.The WRF model overestimates wind speed variability by more than 10% at 925 hPa and 850 hPa, but slightly underestimates variability at higher mandatory levels.
The GEOS-5 analysis underestimates variability consistently at all levels.Both models' correlation skill scores improve with height, as expected (Table 5 and Tables S1-S3).WRF overestimates day-to-day variability, as shown by the variance ratio, in winds at lower levels, and has a larger RMSE than GEOS-5 at lower, but not higher, levels.GEOS-5 underestimates day-to-day variability at all levels.Annual wind direction bias and RMSE, for both models at all sites and mandatory levels are smaller than the number of degrees separating reportable wind directions (not shown).However, at the 925 hPa mandatory level there are seasonal directional biases at many sites in summer and at west coast sites in all seasons (Figure 10).This directional bias may contribute to the summer boundary layer differences in XCO 2 simulations between the two models.Based on these results, the potential improvement from finer resolution WRF simulations in the resolution of mesoscale features and vertical structure near the surface is not evident from our continental-scale evaluations.We discuss further in Section 4 the impact of potential biases in WRF near the surface and in GEOS-5 at upper levels.Based on the spatiotemporal distribution of GOSAT XCO 2 soundings, there is very good agreement between the two modeling systems and GOSAT in the gross seasonal comparisons, well within the uncertainty of the individual GOSAT retrievals (0.5-2.00 ppm, reported with each sounding).Our primary goal of this framework is to control the CO 2 mole fractions introduced into the WRF domain to be as close as possible to the mole fractions from the global model.By scaling the surface fluxes and by constraining the mole fractions in the flow at the boundary walls, we have approached this goal of mass conservation.
The largest spatiotemporal differences in XCO 2 shown in this experiment appear to be due to the model-model differences in transport of the anomalous July biogenic source in the continental northwest and offsetting sinks in the midwest and east in the optimized CMS-Flux biogenic fluxes.However, the very small differences in simulated XCO 2 values mask larger differences within the columns.For example, we can see model-model differences in the transport of the surface flux anomaly in the Pacific Northwest in summer, both within and above the boundary layer (Figure 6).Apparent model differences in vertical mixing within the boundary layer in summer in the eastern part of North America result in lower values within the boundary layer (-0.5 to +0.1 ppm) and higher values above the boundary layer (by +0.1 to +0.5 ppm) in the mesoscale model compared to the global model.We know that this configuration of WRF-Chem lacks convective mass transport of the CO 2 tracers.This will affect the vertical transport of CO 2 into and out of the boundary layer.The CMS-Flux optimzied CO 2 mole fractions also show effects of recent fluxes in the surface layer, and then nearly homogeneous mixing in the next several model layers, as seen in the example profile in Figure 5a.
In this experiment, the meteorological evaluation of the WRF results was approximately comparable to the GEOS-Chem GEOS-5 results, and neither compared well in all respects with the rawinsonde observations.At the selected sites, the GEOS-5 winds were slow by up to 3 m s 1 at nearly all the sites and mandatory levels we used for comparison, while the WRF winds were closer to observations at higher levels.We had hoped to see better performance from the increased resolution, both horizontal and vertical, in WRF.Perhaps the horizontal resolution in WRF is still too coarse to take advantage of any truly mesoscale effects.We also did not assimilate meteorological observations into WRF; this was by design as the primary focus of the experiment was to identify differences in transport of the CO 2 originating from the CMS-Flux optimized biogenic fluxes.
The WRF resolution, the assimilation of meteorological data in WRF, and changes to the boundary layer parameterizations within WRF could be tested to review the conclusions we see here.
One of the main objectives of satellite XCO 2 programs is to provide observations for assimilation into atmospheric inversions to improve the quality of inferred surface fluxes.The satellites provide good coverage of the North American domain in the summer, the season with the most biogenic activity.However, coverage in other seasons is limited, and must be supplemented with other CO 2 observations, forcing the challenge of assimilating both column and surface observations in the same inversion.It is an interesting thought experiment to see what the model-model differences would be if we did have complete satellite sampling coverage over the course of a year.We sampled each model at its own horizontal grid resolution at the hour of the day with the most GOSAT XCO 2 soundings (20 UTC), created simulated XCO 2 from CO 2 columns interpolated to the pressure levels commonly used in the ACOS algorithm, and examined model-model differences in the same way as we did with the simulated GOSAT soundings.We compared the differences for summer, the season with the best GOSAT spatial coverage (See Figures S4-S6).WRF values below 850 hPa are lower than the CMS-Flux values in the east, and WRF values are higher over most of the continent above 850 hPa.These results are reasonably consistent with the summer results shown in Figures 4 and 6, which are conditioned on the spatiotemporal coverage of the GOSAT soundings.This is encouraging.While we will never have perfect coverage in every season with a satellite XCO 2 product, the results presented here show that model-model differences do not appear to be overly dependent on the GOSAT sampling coverage.
We have created a viable framework for comparing the transport of surface fluxes optimized in one model in another model.
Theoretically, this allows for comparisons using different grid resolutions, meteorological drivers, and model parameterization options.In this experiment, we do see some differences between models in horizontal winds and boundary layer mixing, but see no clear advantage of the mesoscale model in the simulation of the satellite-derived XCO 2 .As a caveat, this experiment used one combination of parameterizations of boundary layer processes in the mesoscale model.The framework does provide a computationally feasible laboratory for investigating other possible model configurations.Repeating the experiment in WRF in other configurations (boundary layer physics parameterizations, driver meteorological fields, assimilation of observed meteorology, or finer model resolutions) could yield an ensemble of results that might better establish an envelope of transport uncertainty for the CMS-Flux optimized biogenic flux solution.

Conclusions
We have established a framework for effectively nesting a mesoscale model within a global model achieving approximate mass conservation of the trace gas CO 2 .The CMS-Flux and WRF simulated column-averaged GOSAT XCO 2 samples are comparable within a few parts per million with some spatial differences, across all seasons and all geographic locations in the WRF North American domain.The models used here have very different horizontal and vertical resolutions and computational grids.Scaling the surface fluxes appropriately could be done with a readily-available mass-conserving regridding utility rather than using the custom-computed scaling factors used in this experiment.The code for the provision of the boundary conditions from the global model to the mesoscale model is available from https://github.com/psu-inversion.There is a version for TM5based models such as CarbonTracker (https://www.esrl.noaa.gov/gmd/ccgg/carbontracker),and versions can be created with straightforward modification for other global models with rectangular grids.This framework makes it possible to follow the CO 2 from surface fluxes optimized in a global model into a regional domain, allowing for testing of various transport options.
Hourly WRF 3D CO 2 fields for 2010 from this experiment are archived at the ORNL DAAC (Lauvaux and Butler, 2016) for community use.
Although simulations of entire columns differ very little between the models in this experiment, there are vertical differences in the distribution of the CO 2 within the column, likely attributable to differences in mixing within the two model regimes.
These differences, insignificant when considering total column XCO 2 from the models, will affect the inverse sources and sinks, attributing signals to the surface and the boundaries differently for each model.Unfortunately, both models show significant biases with respect to meteorological observations that vary with season.At the sites we examined, wind speeds are      When released as surface sources, they add  0.03 ppm to the XCO2 reported here.
4. Multiply the mass exchange assigned to each WRF grid cell by the domain-wide scaling factor.These steps are repeated for each of the surface fluxes.The goal is to achieve equal total surface mass exchanges at the hourly resolution of the flux input into the WRF model domain.The component fluxes from CMS-Flux have temporal resolutions Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-342Manuscript under review for journal Geosci.Model Dev. Discussion started: 1 February 2019 c Author(s) 2019.CC BY 4.0 License.varying from hourly to monthly.The optimized biogenic flux is at monthly resolution, but with a 3-hour diurnal cycle overlay with monthly net zero emission/uptake.We do not scale the diurnal cycle overlay.Most monthly scaling factors are in the range [1.0, 1.3] with the exception of some of the minor fluxes (biofuel and ship bunker fuel) with slightly larger scaling factors.The scaled surface fluxes are introduced into the WRF modeling system using the WRFCHEMI function of WRF-Chem.Realistic flux scaling results may also be achieved using a more sophisticated approach, such as the mass-conserving utilities within NCL ESMF (http://www.ncl.ucar.edu/Applications/ESMF.shtml), which would also retain the global model flux patterns, as well as preserving the domain mean.
Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-342Manuscript under review for journal Geosci.Model Dev. Discussion started: 1 February 2019 c Author(s) 2019.CC BY 4.0 License.interpolated optimized mole fractions for both the WRF and CMS XCO 2 computations.For ACOS GOSAT soundings with profile pressure levels below the model surface, we use the CO 2 at the midpoint of the surface model layer.
Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-342Manuscript under review for journal Geosci.Model Dev. Discussion started: 1 February 2019 c Author(s) 2019.CC BY 4.0 License.slightly larger than the mean differences, except for the CMS-Flux simulations in summer and the WRF simulations in summer and fall.Although the model results are centered well with the GOSAT soundings, neither model produces simulations with the full range of values of the GOSAT XCO 2 .We might have expected more variability in the WRF simulations, as multiple Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-342Manuscript under review for journal Geosci.Model Dev. Discussion started: 1 February 2019 c Author(s) 2019.CC BY 4.0 License.

Figure 1 .Figure 2 .Figure 3 .
Figure1.Conceptual illustration of the vertical mapping scheme between an example grid cell from the CMS-Flux GEOS-Chem 47-level grid (left) and the corresponding domain boundary grid cell in the WRF-Chem 50-level grid (right).Numbers indicate the levels in each model.In this case, the mass at the pressure midpoint in level 6 of the WRF column is matched to level 2 in the GEOS-Chem column.With no additional interpolation, level 6 CO2 in the WRF column will be sourced from level 2 in GEOS-Chem (see text).

Figure 4 .Figure 5 .Figure 6 .Figure 8 .
Figure 4. Seasonal mean spatial differences in WRF and CMS simulated XCO2 aggregated to the GEOS-Chem CMS-Flux grid (light gray lines) for (a) spring, (b) summer, (c) fall and (d) winter.Grid cells with no shading have fewer than 10 GOSAT soundings for the season indicated.

*
These non-surface sources are not included in the results presented in this work.

Table 1 .
Surface fluxes for the North American domain in this study

Table 3 .
Model residuals [ppm]relative to ACOSv3.5 GOSAT XCO2.The span between the distribution quantiles Q3 and Q1 defines the interquartile range (IQR).IQRs for residuals of both models and all seasons are in the range [1.8, 1.9] ppm.

Table 5 .
Model comparison to rawinsonde wind speed at mandatory reporting level 925 hPa, 00 UTC in 2010