CO₂ Flux over the Contiguous United States in 2016 Inverted by WRF-Chem/DART from OCO-2 XCO₂ Retrievals

Abstract

High spatial resolution carbon dioxide (CO₂) flux inversion systems are needed to support the global stocktake required by the Paris Agreement and to complement the bottom-up emission inventories. Based on the work of Zhang, a regional CO₂ flux inversion system capable of assimilating the column-averaged dry air mole fractions of CO₂ (XCO₂) retrieved from Orbiting Carbon Observatory-2 (OCO-2) observations had been developed. To evaluate the system, under the constraints of the initial state and boundary conditions extracted from the CarbonTracker 2017 product (CT2017), the annual CO₂ flux over the contiguous United States in 2016 was inverted (1.08 Pg C yr⁻¹) and compared with the corresponding posterior CO₂ fluxes extracted from OCO-2 model intercomparison project (OCO-2 MIP) (mean: 0.76 Pg C yr⁻¹, standard deviation: 0.29 Pg C yr⁻¹, 9 models in total) and CT2017 (1.19 Pg C yr⁻¹). The uncertainty of the inverted CO₂ flux was reduced by 14.71% compared to the prior flux. The annual mean XCO₂ estimated by the inversion system was 403.67 ppm, which was 0.11 ppm smaller than the result (403.78 ppm) simulated by a parallel experiment without assimilating the OCO-2 retrievals and closer to the result of CT2017 (403.29 ppm). Independent CO₂ flux and concentration measurements from towers, aircraft, and Total Carbon Column Observing Network (TCCON) were used to evaluate the results. Mean bias error (MBE) between the inverted CO₂ flux and flux measurements was 0.73 g C m⁻² d⁻¹, was reduced by 22.34% and 28.43% compared to those of the prior flux and CT2017, respectively. MBEs between the CO₂ concentrations estimated by the inversion system and concentration measurements from TCCON, towers, and aircraft were reduced by 52.78%, 96.45%, and 75%, respectively, compared to those of the parallel experiment. The experiment proved that CO₂ emission hotspots indicated by the inverted annual CO₂ flux with a relatively high spatial resolution of 50 km consisted well with the locations of most major metropolitan/urban areas in the contiguous United States, which demonstrated the potential of combing satellite observations with high spatial resolution CO₂ flux inversion system in supporting the global stocktake.

1. Introduction

Since the industrial revolution, the significant increase of carbon dioxide (CO₂) in the atmosphere caused by anthropogenic activities is believed to be one of the main driving forces of climate change [1]. The atmospheric CO₂ concentration has increased from 277 parts per million (ppm) in 1750 to above 409 ppm in 2019 [2].

In order to mitigate climate change, the Paris Agreement [3], a legally binding international treaty on climate change, was adopted at the 21st Conference of the Parties of the United Nations Framework Convention on Climate Change (UNFCCC) and entered into force on November 4, 2016. The Paris Agreement [3] requires all parties to submit their emission reduction plans in the form of nationally determined contributions (NDC) and requires regular global stocktake processes to verify whether all parties have fulfilled their commitments.

Currently, the statistical bottom-up methodology developed by the Inter-governmental Panel on Climate Change (IPCC) is used to prepare the emission inventory reports required by UNFCCC [4]. This methodology uses statistical activity data and source-specific emission factors to quantify the anthropogenic CO₂ emission of a country or a region. However, due to large uncertainties in the emission factors and in the statistical activity data, as well as biases caused by unaccounted sources, the bottom-up emission inventories can have significant uncertainties, especially for countries with underdeveloped statistical infrastructure. Therefore, with the progress in observations and inverse models, the Subsidiary Body for Scientific and Technological Advice (SBSTA) of UNFCCC [5,6] along with the IPCC Task Force [4] acknowledged the complementary capability offered by the top-down observation approach, which links emissions with atmospheric concentration observations by atmospheric transport and chemistry models, i.e., the inverse models.

Before satellite-based observations of atmospheric CO₂ concentrations were widely available, the research of top-down approach mainly focused on inversing CO₂ flux from ground-based observations [7,8,9,10,11,12]. Compared with ground-based observations, satellite observations have the advantage of providing much better spatiotemporal coverage. With the development of technology, in the past 10 years or so, several CO₂ monitoring satellites—including the Greenhouse gases Observing SATellite (GOSAT) [13], GOSAT-2 [14], the Orbiting Carbon Observatory 2 (OCO-2) [15], OCO-3 [16,17], and TanSat [18], had been launched. Several others, such as MicroCarb [19,20], geoCARB [21,22], and a CO2M constellation [23,24]—are currently under development.

Promoted by the CO₂ monitoring satellites, recent studies used the column-averaged dry air mole fractions of CO₂ (XCO₂) products retrieved from satellite observations to inverse CO₂ flux [25,26,27,28,29,30]. Among them, Deng et al. [31] used the four-dimensional variational (4D-Var) approach coupled with GEOS-Chem model to demonstrate the significant impact of GOSAT observations on the flux inversion, Crowell et al. [32] utilized different atmospheric transport models, data assimilation methods and prior fluxes to inverse the global CO₂ fluxes from OCO-2 observations and evaluate their uncertainties, Chevallier, et al. [33] and Wang et al. [34] assimilated observations from GOSAT and OCO-2 respectively to conduct global CO₂ flux inversion. Although these global-scale inversion studies can deepen the understanding of the global carbon cycle mechanism that affects climate change, due to their coarse spatiotemporal resolution, their results are not enough to support the global stocktake required by the Paris Agreement.

To complement the bottom-up emission inventories, regional flux inversion systems with high spatiotemporal resolution are needed. Villalobos et al. [35] used the Community Multiscale Air Quality (CMAQ) transport-dispersion model coupled with 4D-Var to study the impact of OCO-2 observations on CO₂ flux inversions and on their uncertainty reduction over Australia through a series of Observing System Simulation Experiments (OSSEs). Zheng et al. [36] and Peng et al. [37] developed regional CO₂ flux inversion systems based on WRF and CMAQ modeling system coupled with 4D-Var and ensemble Kalman filter approach, and then evaluated the system performances in the continental United States and East Asia through OSSEs, respectively. However, real XCO₂ retrievals from OCO-2 were not yet utilized in these studies.

In this study, based on the work of Zhang et al. [38], a regional CO₂ flux inversion system had been developed from WRF-Chem model [39] coupled with Data Assimilation Research Testbed (DART) [40] using the ensemble adjustment Kalman filter (EAKF) [41,42] assimilation method to invert regional CO₂ flux from OCO-2 XCO₂ retrievals. The contiguous United States was selected as the study area because there were relatively abundant ground-based observations available for evaluation. The study period was the entire year of 2016. Different types of CO₂ concentration observations, CO₂ flux measurements, as well as the products from the OCO-2 model intercomparison project (MIP) [32] and CT2017 [10], were used to evaluate the performance of the inversion system over the contiguous United States.

2. Materials and Methods

2.1. CO₂ Transport Model

WRF-Chem model version 3.9.1 [39] was used as the CO₂ transport model. As shown in Figure 1, the model domain, having 103 (west–east) × 82 (south–north) grid points with a spatial resolution of 50 km, covered the contiguous United States. In the vertical direction, the atmosphere was decomposed into 45 terrain-following levels from the surface to the upper boundary at 50 hPa. The physical and chemical options of the model, listed in

Table 1. WRF-Chem configurations.

Model Options	Configurations
Domain center	34.939 °N, −96.275 °W
Grid resolution	50 km
nx, ny, nz	103, 82, 45
Time step	240 s
Microphysics process	WSM 5-class simple ice scheme
Cumulus parameterization	Kain-Fritsch scheme
Longwave atmospheric radiation	RRTM scheme
Shortwave atmospheric radiation	Dudhia scheme
Planetary boundary layer scheme	MYNN 2.5 level TKE
Surface layer scheme	MYNN
Land surface scheme	Unified Noah Land surface model
Chemical option	chem_opt = 16 (CO2 only)

, were the same as those used in Zhang, et al. [38]. The chemical option “chem_opt = 16” means that CO₂ is treated as a passive tracer with no feedback on meteorological variables and is not part of any chemical reactions.

The Final (FNL) Operational Global Analysis data [43] with 1° × 1° spatial resolution and 6-h interval from National Centers for Environmental Prediction (NCEP), were used as the meteorological initial and boundary conditions of the model. The initial and boundary conditions of CO₂ concentrations were interpolated from the CO₂ mole fraction product of the NOAA’s CarbonTracker, version 2017 (CT2017) [10]. Prior CO₂ fluxes used by CarbonTracker in CT2017, including anthropogenic emissions, fire emissions, biogenic fluxes, and ocean fluxes were used as the prior fluxes input to the model.

2.2. OCO-2 XCO₂ Retrievals

OCO-2 is an atmospheric CO₂ observing satellite launched by the National Aeronautics and Space Administration (NASA) in 2014 [44,45]. OCO-2 is routinely making around 1 million observations each day, and more than 10% of these observations are sufficiently cloud-free to allow for the retrieval of XCO₂ via an inversion algorithm [46]. The XCO₂ retrieved from OCO-2 observations have a median difference of less than 0.5 ppm and root mean square differences typically below 1.5 ppm compared with Total Carbon Column Observing Network (TCCON) [47].

The Version 9r of Level 2 Lite products of the OCO-2 XCO₂ retrievals [48,49] were used in this study. In order to be consistent with the OCO-2 MIP [32], in the lite products, only the retrievals from nadir observation mode with good quantity based on the “xco2_quality_flag” were selected. Rather than assimilating each XCO₂ retrieval falling inside a model grid cell separately, a representative mean XCO₂ retrieval and its associated uncertainty of a model grid cell were calculated and assimilated according to the strategy of Crowell et al. [32].

The OCO-2 XCO₂ retrievals over the contiguous United States in each month of 2016 were processed by this strategy and the results are shown in Figure 2.

2.3. Observation Operator

The CO₂ concentrations forecasted by WRF-Chem in the model space were converted into XCO₂ values in the OCO-2 observation space by an observation operator adopted from Connor et al. [50]. The observation operator $H$was defined as

$$XC{O}_2^m=XC{O}_2^b+{\displaystyle \sum }_j{h}_j{a}_{C{O}_2,j}{\left(x^m-x^b\right)}_j$$

(1)

where $XC{O}_2^b$, $h$, $a_{C{O}_2}$, and $x^b$ are the prior XCO₂ value, the pressure weighting function, the column averaging kernel, and the prior CO₂ concentration profile used by OCO-2 XCO₂ retrieval processes, respectively. All these parameters are provided in the OCO-2 Level 2 Lite products and require pretreatments similar to the representative mean XCO₂ retrievals. $x^m$ is the vertical profile of CO₂ concentrations interpolated from WRF-Chem forecasts to the pressure levels of the OCO-2 XCO₂ retrieval. $XC{O}_2^m$ is the converted result of the WRF-Chem forecasts by the observation operator $H$.

2.4. Regional CO₂ Flux Inversion System

Zhang et al. [38] had implemented the assimilation of the OCO-2 XCO₂ retrievals by extending DART [40] to improve the estimation of regional CO₂ concentrations. DART [40] is a widely used open-source software framework for ensemble data assimilation research [51,52,53,54].

On the basis of Zhang’s assimilation system [38], a regional flux inversion system had been further developed in this study by augmenting the state vector to include the CO₂ fluxes [55,56]. The state vector $x$ was defined as

$$x={\left[CC,CF\right]}^T$$

(2)

where $CC$ represent the CO₂ concentration states in the model space, which was advanced by WRF-Chem and $CF$represent the CO₂ fluxes input into WRF-Chem. $CF$ was the non-observable part of the state vector and could not be forwarded by WRF-Chem model.

The initial value of $CC$ for the $k$-th ensemble member, defined as $C{C}_k^{init}$, was generated in the same way as that of Zhang, et al. [38] by imposing zero-mean Gaussian distributions to the average CO₂ concentrations interpolated from CT2017 [10]. The standard deviations of these Gaussian distributions were assumed to be 1.5% of the average CO₂ concentrations.

The update process of $CC$ was the same as that of Zhang et al. [38]. Let $C{C}_{k,i}^f$ and $C{C}_{k,i}^a$ represent the forecast of $CC$ advanced by WRF-Chem starting from $C{C}_k^{init}$ and the analysis result of $CC$ for the $k$-th ensemble member at time $i \left(i\ge 1\right)$, respectively. The $C{C}_{k,i}^f$ was converted into the observation space by the observation operator $H$ defined by Equation (3), i.e., $y_{k,i}^f=H\left(C{C}_{k,i}^f\right)$. Then the analysis result $y_{k,i}^a$ of the $k$-th ensemble member at time $i$ was calculated by

$$y_{k,i}^a=\left[\sqrt{\frac{{\widetilde{{\sigma }_{XC{O}_2}}}^2 }{{\widetilde{{\sigma }_{XC{O}_2}}}^2+{\overline{\sigma }}_i^f{}^2}}\right]\left[y_{k,i}^f-{\overline{y}}_i^f\right]+\left[\frac{{\overline{y}}_i^f}{{\overline{\sigma }}_i^f{}^2}+\frac{\widetilde{XC{O}_2}}{{\widetilde{{\sigma }_{XC{O}_2}}}^2 }\right]{\left[\frac{1}{{\overline{\sigma }}_i^f{}^2}+\frac{1}{{\widetilde{{\sigma }_{XC{O}_2}}}^2 }\right]}^{-1}$$

(3)

where ${\overline{y}}_i^f$ and ${\overline{\sigma }}_i^f$ were the forecast mean and the spread of all ensemble members in the observation space at time $i \left(i\ge 1\right)$, respectively. $\widetilde{XC{O}_2}$ and $\widetilde{{\sigma }_{XC{O}_2}}$ were the representative mean XCO₂ retrieval and its associated uncertainty of a model grid cell calculated from the OCO-2 retrievals according to the strategy of Crowell et al. [32].

The analysis result of $C{C}_{k,i}^f$ was updated as

$$C{C}_{k,i}^a=C{C}_{k,i}^f+\alpha \frac{\sigma \left(C{C}_i^f,y_i^f\right)}{{\overline{\sigma }}_i^f{}^2}\left(y_{k,i}^a-y_{k,i}^f\right)$$

(4)

where $\alpha$ was a covariance localization coefficient to compensate the sampling error caused by the limited ensemble size [57], $\sigma \left(C{C}_i^f,y_i^f\right)$ is the covariance of the forecasts of $CC$ and their corresponding conversion results in the observation space across all ensemble members at time $i$ $\left(i\ge 1\right)$.

The initial value of $CF$ for the $k$-th ensemble member at time $i$ $\left(i\ge 1\right)$, defined as $C{F}_{k,i-1}^{init}$, was generated by imposing zero-mean Gaussian distributions to the average CO₂ fluxes interpolated from the prior CO₂ fluxes used by CT2017 [10]. The standard deviations of these flux Gaussian distributions were assumed to be 50% of the average CO₂ fluxes similar to Chevallier et al. [30].

Let $C{F}_{k,i-1}^{prior}$and $C{F}_{k,i-1}^{post}$ represent the prior and posterior CO₂ fluxes of the model grids for the $k$-th ensemble member at time $\left(i\ge 1\right)$, respectively. Let ${\lambda }_{k,i-1}^{prior}$ and ${\lambda }_{k,i-1}^{post}$ represent the prior and posterior scaling factors between the pairs of $C{F}_{k,i-1}^{prior}$, $C{F}_{k,i-1}^{init}$ and $C{F}_{k,i-1}^{post}$, $C{F}_{k,i-1}^{prior}$, respectively. The $C{F}_{k,i-1}^{prior}$ was calculated as

$$C{F}_{k,i-1}^{prior}={\lambda }_{k,i-1}^{prior}\times C{F}_{k,i-1}^{init} ,i\ge 1$$

(5)

where the prior scaling factor ${\lambda }_{k,i-1}^{prior}$ was evolved forward with time $i$ by

$${\lambda }_{k,i-1}^{prior}=\left\{\begin{array}{c} 1 , i=1\\ {\lambda }_{k,i-2}^{post} , i\ge 2\end{array}\right.$$

(6)

The posterior CO₂ flux and scaling factor were updated by

$$C{F}_{k,i-1}^{post}=C{F}_{k,i-1}^{prior}+\alpha \frac{\sigma \left(C{F}_{i-1}^{prior},y_i^f\right)}{{\overline{\sigma }}_i^f{}^2}\left(y_{k,i}^a-y_{k,i}^f\right) ,i\ge 1$$

(7)

$${\lambda }_{k,i-1}^{post}=C{F}_{k,i-1}^{post}/C{F}_{k,i-1}^{prior},i\ge 1$$

(8)

where $\alpha$ and $\sigma \left(C{F}_{i-1}^{prior},y_i^f\right)$ have the similar definitions as those of in the Equation (4).

2.5. Experiment Design

Two experiments were conducted to evaluate the performance of the flux inversion system. The first one, named “DA_FLUX”, used the flux inversion system to optimize both the atmospheric CO₂ concentrations and CO₂ fluxes during the assimilation process. The second experiment, named “SIM”, only used WRF-Chem to simulate the atmospheric CO₂ concentrations and did not assimilate any OCO-2 XCO₂ retrievals.

2.5.1. The DA_FLUX Experiment

In the DA_FLUX experiment, the CO₂ concentrations and CO₂ fluxes over the contiguous United States in 2016 were estimated by the flux inversion system through the assimilation of OCO-2 XCO₂ retrievals.

The ensemble size of the system was 20. The initial and boundary conditions of CO₂ concentrations were interpolated from CT2017 [10]. A strategy similar to CarbonTracker was taken to quantify the uncertainty due to prior flux selection (https://www.esrl.noaa.gov/gmd/ccgg/carbontracker/CT2019B_doc.php, last access: 18 June 2021), by interpolating two sets of prior fluxes from two different combinations of terrestrial fluxes, wildfire emissions, air–sea CO₂ exchanges, and fossil fuel emissions from CT2017. The mean posterior CO₂ fluxes of two inversions, each using one set of prior fluxes, were the final CO₂ flux results of the inversion system. The final CO₂ concentration results were estimated by applying the mean posterior fluxes to the forward model WRF-Chem once again and without assimilation.

The study period was the entire year of 2016. The period of five weeks before 1 January 2016 was used as the spin-up phase of the simulation. The simulation cycle was 6-h at 00:00, 06:00, 12:00, and 18:00 UTC, respectively. A short assimilation window of one day was used according to Liu et al. [58].

2.5.2. The SIM Experiment

The SIM experiment was a control experiment, using only WRF-Chem to simulate the CO₂ concentrations, and not assimilating any OCO-2 XCO₂ retrievals. The prior CO₂ flux of the SIM experiment was the average of the two sets of prior fluxes used in the DA_FLUX experiment. All other initial conditions and configurations were the same as those of the DA_FLUX experiment.

2.6. Evaluation

The level-4 flux product version 7 (https://www.esrl.noaa.gov/gmd/ccgg/OCO2/, last access: 18 June 2021) from the OCO-2 model intercomparison project (OCO-2 MIP) was used to compare with the posterior CO₂ flux obtained by the regional flux inversion system. The OCO-2 MIP [32] is a collaboration organized by the OCO-2 Science Team to study the impact of assimilating OCO-2 XCO₂ retrievals into atmospheric inversion models. Currently, there are nine models in the OCO-2 MIP version 7, they differ in the sources of prior CO₂ flux, the types of atmospheric chemical transport model, the methods of data assimilation, and so on. The detailed information of these models is listed in Table S1.

Independent ground and aircraft observations were used to evaluate the performance of the regional flux inversion system. The CO₂ flux observations from AmeriFlux [59] and ONEFlux [60], including 18 and 10 eddy covariance flux measurement sites respectively, were used to evaluate the CO₂ flux results of the DA_FLUX experiment. The detailed information of these flux observation sites is listed in Tables S2 and S3, and their locations are shown in Figure 1e,f, respectively.

The XCO₂ retrievals of four TCCON [61] sites, including Park Falls [62], Lamont [63], Edwards [64], and Caltech [65], CO₂ concentration observations of 10 tower sites and 20 aircraft campaigns from the GLOBALVIEWplus project [66] with Observation Package (ObsPack) [67] data products “obspack_co2_1_GLOBALVIEWplus_v5.0_2019-08-12” were used to evaluate the CO₂ concentration results of the two experiments. The detailed information of these sites and campaigns is listed in Tables S4 and S5, and their locations are shown in Figure 1a,d.

The posterior CO₂ concentrations and posterior fluxes from CT2017 were also compared with the independent observations to serve as references for the evaluation.

The experimental results were evaluated by mean bias error (MBE), root-mean-square error (RMSE), and correlation coefficient (CORR), same as the Equations (6), (8), and (9) in Zhang’s work [38].

The uncertainty reduction between the posterior and prior CO₂ fluxes is a useful method to evaluate the impact of assimilation of OCO-2 XCO₂ retrievals in constraining CO₂ flux. The uncertainty reduction achieved by the flux inversion system was evaluated by the equation [68]

$$UR=\left(1-\frac{{\sigma }_{post}}{{\sigma }_{prior}}\right)\times 100\%$$

(9)

where ${\sigma }_{post}$ and ${\sigma }_{prior}$ were the standard deviations of the posterior and prior CO₂ fluxes, respectively

3. Results and Discussion

3.1. Reginal CO₂ Flux Inversion Results

By assimilating the OCO-2 XCO₂ retrievals, the regional CO₂ flux inversion system in the DA_FLUX experiment optimized the annual mean CO₂ flux over the contiguous United States in 2016 from the prior flux of 1.51 ± 0.11 Pg C yr⁻¹ to the posterior flux of 1.08 ± 0.03 Pg C yr⁻¹, which was more in line with the fluxes estimated by CT2017 and the models from OCO-2 MIP. The prior and posterior monthly mean CO₂ fluxes of the DA_FLUX experiment are listed in

Table 2. Prior and posterior monthly mean CO₂ fluxes of the DA_FLUX experiment over the contiguous United States in 2016.

Month	Prior Fluxes (Pg C)	Posterior Fluxes (Pg C)
1	0.35 ± 0.03	0.33 ± 0.01
2	0.32 ± 0.008	0.27 ± 0.001
3	0.28 ± 0.02	0.19 ± 0.02
4	0.12 ± 0.05	0.05 ± 0.03
5	−0.05 ± 0.06	−0.1 ± 0.05
6	−0.21 ± 0.05	−0.17 ± 0.05
7	−0.20 ± 0.05	−0.07 ± 0.02
8	−0.08 ± 0.002	−0.02 ± 0.001
9	0.04 ± 0.03	0.07 ± 0.02
10	0.21 ± 0.03	0.16 ± 0.03
11	0.35 ± 0.01	0.19 ± 0.03
12	0.37 ± 0.02	0.19 ± 0.04
Annual	1.51 ± 0.11	1.08 ± 0.03

. The results of DA_FLUX indicated that in 2016 the contiguous United States was a carbon sink during the four-month period from May to August due to the vegetation growing season, and was a carbon source in the remaining eight months.

The spatial distribution of the annual mean posterior CO₂ flux of DA_FLUX, shown in Figure 3, indicates the areas with higher CO₂ fluxes roughly coincided with the densely populated urban areas—i.e., the metropolis areas of New York, New Jersey, Detroit of Michigan, Chicago of Illinois, Atlanta of Georgia, Tampa and Orlando of Florida, Miami of Florida, Houston of Texas, Los Angeles, and San Francisco of California all showed up as CO₂ sources of varying magnitude. It implies the majority of the CO₂ emission in the contiguous United States comes from the entire northeast, midwest, and a few hotspots in the south and west because these are the largest metropolitan/urban areas. Figure 3 also indicates almost the entire Great Plains and Rocky Mountain West were a weak CO₂ sink because with a few exceptions there are virtually no large cities. The spatial distributions of the monthly mean posterior fluxes are shown in Figure S1.

3.1.1. Compared with CT2017 and OCO-2 MIP Models

Figure 4 and Figure 5 show the comparison of the annual and monthly posterior CO₂ flux of the DA_FLUX experiment with the CO₂ fluxes estimated by CT2017 and OCO-2 MIP models, and with the prior flux used in the experiment, respectively. The detailed information of the OCO-2 MIP models shown in Figure 4 and Figure 5 is listed in Table S1. The monthly mean values of each flux in Figure 5 are listed in Table S6.

The flux of 1.19 Pg C yr⁻¹ of CT2017 was extracted from the posterior CO₂ flux of CT2017 in 2016 according to the border of the contiguous United States. The fluxes of OCO-2 MIP models were extracted from the “TransCom 02” region of the OCO-2 MIP level 4 flux products, according to the area ratio of the contiguous United States and the “TransCom 02” region. The “TransCom 02” region covers the North American temperate zone, which is slightly larger than the contiguous United States. The mean flux and standard deviation of all the nine models from OCO-2 MIP were 0.76 Pg C yr⁻¹ and 0.29 Pg C yr⁻¹, respectively.

Figure 4 shows the annual posterior CO₂ flux of the DA_FLUX experiment was 0.11 Pg C yr⁻¹ lower than that of CT2017, and was 0.05 Pg C yr⁻¹ higher than the highest flux from OCO-2 MIP estimated by the model of Baker_PCTM_Ldwo_LN and Baker_PCTM_Taka_LN.

Figure 5 shows that the monthly mean fluxes of DA_FLUX from May to July, when the contiguous United States was a carbon sink, were significantly smaller than that of OCO-2 MIP models, i.e., less carbon uptake. This explains why the annual flux of DA_FLUX was higher than that of OCO-2 MIP models. By analyzing the prior flux, the CO₂ concentrations simulated by WRF-Chem before assimilation and the value of representative mean XCO₂ retrievals of OCO-2, it was found that the small posterior fluxes in May and June were mainly caused by the significantly smaller prior fluxes which implied less CO₂ uptake compared with those of OCO-2 MIP models. The prior fluxes in May and June were only 24.3% and 54.8% of the mean fluxes of OCO-2 MIP models, respectively. The small posterior flux in July was caused by the combination of the relatively high representative mean XCO₂ retrievals of OCO-2 and a small prior flux. The average value of the representative mean XCO₂ retrievals in July was 401.21 ppm and was 0.35 ppm higher than the monthly mean CO₂ concentrations simulated by WRF-Chem before assimilation. The prior flux in July was 66.5% of the mean fluxes of OCO-2 MIP models. Figure 5 also shows CT2017 and OCO-2 MIP models were inconsistent in terms of whether the contiguous United States were the sink or source of CO₂ in the two months of April and August. On this issue, the results of the DA_FLUX experiment were consistent with the majority.

Based on the above annual and monthly comparison results with other models, the CO₂ flux estimated by the inversion system presented here can be considered reasonable.

3.1.2. Uncertainty Reduction by the Inversion System

The uncertainty of the inverted annual CO₂ flux in the DA_FLUX experiment was reduced by 14.71% compared to the prior flux.

Figure 2 and Figure 6 show the spatial distribution of the representative mean XCO₂ retrievals of OCO-2 and the mean uncertainty of the prior CO₂ flux of each month in 2016 used in the DA_FLUX experiment, respectively. The monthly mean uncertainty of the prior CO₂ flux, the monthly mean uncertainty reduction percentage over all grid cells, and the number of the representative mean XCO₂ retrievals of OCO-2 of each month in 2016 are listed in Error! Reference source not found. Figure 7 shows the spatial distribution of monthly mean uncertainty reduction between the posterior and the prior CO₂ flux in the DA_FLUX experiment. In Figure 7, areas with negative uncertainty reduction were masked as zero, since such uncertainty increases simply resulted from the small number of realizations [35].

From

Table 3. Monthly mean uncertainty of the prior CO₂ flux, the monthly mean uncertainty reduction percentage over all grid cells and the number of the representative mean XCO₂ retrievals of OCO-2 of each month in 2016.

Month in 2016	Percentage of Mean Uncertainty Reduction over All Grid Cells	Number of Representative Mean XCO2 Retrievals of OCO-2	Mean Uncertainty of the Prior CO2 Flux over All Grid Cells (g C m−2 d−1)
1	14.82%	303	0.16
2	12.91%	458	0.16
3	12.80%	347	0.17
4	19.59%	449	0.22
5	18.30%	504	0.22
6	36.17%	695	0.31
7	38.42%	727	0.32
8	30.18%	429	0.22
9	24.42%	620	0.21
10	21.72%	575	0.17
11	10.33%	566	0.14
12	9.53%	344	0.14
Annual	14.71%	6017	0.085

, the uncertainty reduction had a strong positive correlation with the uncertainty of the prior CO₂ flux and the number of the representative mean XCO₂ retrievals of OCO-2. The correlations between the uncertainty reduction and the prior flux uncertainty, the number of the representative mean XCO₂ retrievals were 0.93 and 0.72, respectively. The largest uncertainty reductions of 36.17% and 38.42% occurred in June and July, which coincided with the period that had the largest prior flux uncertainties of 0.31 g C m⁻² d⁻¹ and 0.32 g C m⁻² d⁻¹, respectively. June and July are the middle of the growing season in the contiguous United States, and their mean prior flux uncertainties over all the grid cells are driven up by vegetation photosynthesis.

In the spatial dimension, except for September, the regions with larger uncertainty reduction in Figure 7 roughly coincided with the locations with higher prior uncertainty in Figure 6. In August and September, Error! Reference source not found. shows that the areas, mainly near the southwest of the Great Lakes of North America, had similar distributions of high prior uncertainty. This “upper midwest” region of the contiguous United States is also referred to as the “corn belt”, it is a known mega-sink of CO₂ during the agricultural growing season [69]. Figure 7 shows this region had a larger uncertainty reduction in August, which was in line with the expectation, while it only had a limited uncertainty reduction in September. By comparing the distribution of OCO-2 representative XCO₂ retrievals in August and September shown in Figure 2, it was found that there were no OCO-2 XCO₂ retrievals in this “upper midwest” region in September. This explains why the region had a limited uncertainty reduction in September, and also demonstrated the impact of OCO-2 XCO₂ retrievals on the flux constraints.

3.1.3. Compared with CO₂ Flux Measurements

Although the spatial representativeness of the CO₂ flux tower measurements is usually less than 1 km², and the grid cell size in this study was 50 km, meaning that there was some mismatch of the spatial scale, the statistical results of comparing the CO₂ flux measurements with the DA_FLUX experimental results were still meaningful.

The annual mean of the prior, posterior CO₂ flux of DA_FLUX and the posterior flux of CT2017 interpolated to the locations of flux measurement sites are listed in

Table 4. Annual mean value of the prior (PRIOR), posterior (DA_FLUX) CO₂ flux of the DA_FLUX experiment, the posterior CO₂ flux of CT2017 interpolated to the locations of flux measurement sites, and annual mean CO₂ flux measurements of AmeriFlux and ONEFlux.

Dataset	Flux Measurements (g C m−2 d−1)	Flux
		PRIOR (g C m−2 d−1)	CT2017 (g C m−2 d−1)	DA_FLUX (g C m−2 d−1)
AmeriFlux	−0.25	0.71	0.85	0.46
ONEFlux	−0.52	0.36	0.35	0.24
All	−0.38	0.54	0.60	0.35

, and box-plotted on the right side of Figure 8. The overall annual mean flux measurement of AmeriFlux [59] and ONEFlux [60] was −0.38 g C m⁻² d⁻¹. The corresponding annual mean posterior flux of DA_FLUX was 0.35 g C m⁻² d⁻¹, which were 0.19 g C m⁻² d⁻¹ and 0.25 g C m⁻² d⁻¹ smaller than that of prior flux and CT2017, respectively, and closer to the measurements.

The line charts on the left side of Figure 8 show the monthly mean of the prior and posterior CO₂ fluxes of DA_FLUX, of the posterior CO₂ flux of CT2017, and of the CO₂ flux measurements. The detailed monthly mean CO₂ flux values are listed in Table S7.

The statistical evaluation results of comparing the CO₂ flux measurements from AmeriFlux [59] and ONEFlux [60] with the prior, posterior CO₂ fluxes of the DA_FLUX experiment, and with the posterior CO₂ flux of CT2017 are listed in

Table 5. Statistical results of comparing the CO₂ flux measurements from AmeriFlux and ONEFlux with the prior (PRIOR), posterior (DA_FLUX) CO₂ flux of the DA_FLUX experiment, and with the posterior CO₂ flux of CT2017.

Dataset	Flux	RMSE (g C m−2 d−1)	MBE (g C m−2 d−1)	CORR
AmeriFlux	PRIOR	0.98	0.97	0.97
	CT2017	1.15	1.10	0.92
	DA_FLUX	0.80	0.71	0.95
ONEFlux	PRIOR	0.94	0.88	0.85
	CT2017	0.94	0.86	0.75
	DA_FLUX	0.83	0.76	0.83
All dataset	PRIOR	0.95	0.94	0.97
	CT2017	1.06	1.02	0.92
	DA_FLUX	0.79	0.73	0.94

. The RMSE and MBE of the posterior CO₂ flux in DA_FLUX against all the flux measurements from AmeriFlux and ONEFlux were 0.79 g C m⁻² d⁻¹ and 0.73 g C m⁻² d⁻¹, respectively, which were 0.16 g C m⁻² d⁻¹ (16.84%) and 0.21 g C m⁻² d⁻¹ (22.34%) smaller than that of the prior flux and were 0.27 g C m⁻² d⁻¹ (25.47%) and 0.29 g C m⁻² d⁻¹ (28.43%) smaller than that of CT2017, respectively.

The evaluation results indicated that the posterior CO₂ fluxes of the inversion system were significantly more consistent with the flux measurements than the prior CO₂ fluxes and the posterior CO₂ flux of CT2017.

3.2. CO₂ Concentrations Results

The annual mean XCO₂ over the contiguous United States in 2016 estimated by the SIM, DA_FLUX experiment, and CT2017 were 403.78 ppm, 403.67 ppm, and 403.29 ppm, respectively, and their distributions are shown in Figure 9. The annual mean XCO₂ of DA_FLUX was 0.11 ppm smaller than that of SIM and closer to CT2017. The differences of annual mean XCO₂ between SIM and CT2017, between DA_FLUX and CT2017, between DA_FLUX and SIM were 0.49 ppm, 0.38 ppm, −0.11 ppm, respectively. Their distributions are shown in Figure 10. It is obvious that the result of DA_FLUX was more in line with CT2017 than that of SIM.

The boxplots of the monthly mean CO₂ concentrations of the SIM, DA_FLUX experiment and CT2017 compared with monthly mean observations of TCCON, tower, and aircraft are shown in Figure 11, the corresponding data of Figure 11 are listed in Table S8. The detailed evaluation results of statistical metrics are shown in

Table 6. Evaluation results between the monthly mean CO₂ concentration or XCO₂ observations and the corresponding results from SIM, DA_FLUX and CT2017, where the CORRs are statistically significant at the level of p < 0.01.

Observation Type	Experiment	RMSE (ppm)	MBE (ppm)	CORR
TCCON	SIM	0.70	0.36	0.95
	DA_FLUX	0.81	0.17	0.92
	CT2017	0.69	−0.11	0.94
Tower	SIM	2.83	1.41	0.92
	DA_FLUX	2.63	−0.05	0.89
	CT2017	2.71	0.70	0.89
Aircraft	SIM	0.95	−0.20	0.94
	DA_FLUX	0.68	−0.05	0.97
	CT2017	0.65	−0.10	0.97

. The comparison of the monthly mean CO₂ concentrations of the SIM, DA_FLUX experiment, and CT2017 with monthly mean observations of TCCON, tower, and aircraft are shown in Figures S2–S4, respectively.

The RMSE between the monthly mean CO₂ concentrations of DA_FLUX and the observations of tower and aircraft were 2.63 ppm and 0.68 ppm, which were 0.20 ppm and 0.27 ppm smaller than those of SIM, respectively. The RMSE between the monthly mean XCO₂ of DA_FLUX and the TCCON observations was 0.81 ppm, which was 0.11 ppm larger than that of SIM. The reason is, shown in Figure S2, the monthly mean XCO₂ values estimated by DA_FLUX spread more widely but symmetrically around the TCCON observations than those of SIM. The MBE between the monthly mean CO₂ concentrations of DA_FLUX and the observations of TCCON, tower, and aircraft were 0.17 ppm, −0.05 ppm, and −0.05 ppm, which were 0.19 ppm, 1.36 ppm, and 0.15 ppm smaller than those of SIM, respectively. Figure S3 and S4 show the comparison of the monthly mean tower and aircraft observations with the modeled CO₂ concentrations from the SIM, DA_FLUX experiment, and CT2017, respectively.

4. Conclusions

A regional CO₂ flux inversion system had been developed from WRF-Chem model coupled with DART using the EAKF assimilation method to invert regional CO₂ flux from OCO-2 XCO₂ retrievals. The effectiveness of the inversion system was verified by the DA_FLUX experiment over the contiguous United States in 2016. The posterior CO₂ fluxes of the DA_FLUX experiment were compared with the corresponding fluxes extracted from OCO-2 MIP and CT2017 and was evaluated with the independent ground-based flux measurements. The inverted annual CO₂ flux over the contiguous United States in 2016 was 1.08 Pg C yr⁻¹, which was 0.05 Pg C yr⁻¹ higher than the largest flux from two models in OCO-2 MIP and was 0.11 Pg C yr⁻¹ lower than the posterior CO₂ flux of CT2017. The uncertainty of the inverted annual CO₂ flux was reduced by 14.71% compared to the prior flux, which demonstrated the impact of OCO-2 XCO₂ retrievals on the flux constraints. The evaluation with the flux measurements showed that RMSE and MBE between the posterior CO₂ flux of the DA_FLUX experiment and the ground observations were reduced by 16.84% and 22.34% compared with that of the prior flux and were reduced by 25.47% and 28.43% compared with that of CT2017, respectively. The CO₂ concentrations estimated by the DA_FLUX experiment were more in line with CT2017 than that of the SIM experiment which did not assimilate the OCO-2 XCO₂ retrievals. The annual mean XCO₂ of the DA_FLUX experiment over the contiguous United States in 2016 was 403.67 ppm which was 0.11 ppm smaller than that of the SIM experiment (403.78 ppm) and closer to CT2017 (403.29 ppm). MBE between the inverted CO₂ flux and flux measurements was 0.73 g C m⁻² d⁻¹, was reduced by 22.34% and 28.43% compared to those of the prior flux and CT2017, respectively. MBEs between the CO₂ concentrations of the DA_FLUX experiment and concentration measurements from TCCON, towers, and aircraft were reduced by 52.78%, 96.45%, and 75%, respectively, compared to those of the SIM experiment. These experiments and evaluations indicated that the regional CO₂ flux inversion system had a reasonable performance.

Finally, the significant drop in uncertainty reduction of the inverted CO₂ flux of the “Upper Midwest” region in September 2016 due to the lack of OCO-2 XCO₂ retrievals, and the fact that CO₂ emission hotspots indicated by the inverted annual CO₂ flux with a relatively high spatial resolution of 50 km consisted well with the locations of most major metropolitan/urban areas in the contiguous United States, demonstrated the indispensable impact of OCO-2 XCO₂ retrievals on the regional flux constraints and the potential of combining satellite observations with high spatial resolution CO₂ flux inversion system in supporting the global stocktake.