Regional CO2 Inversion Through Ensemble‐Based Simultaneous State and Parameter Estimation: TRACE Framework and Controlled Experiments

Atmospheric inversions provide estimates of carbon dioxide (CO2) fluxes between the surface and atmosphere based on atmospheric CO2 concentration observations. The number of CO2 observations is projected to increase severalfold in the next decades from expanding in situ networks and next‐generation CO2‐observing satellites, providing both an opportunity and a challenge for inversions. This study introduces the TRACE Regional Atmosphere–Carbon Ensemble (TRACE) system, which employ an ensemble‐based simultaneous state and parameter estimation (ESSPE) approach to enable the assimilation of large volumes of observations for constraining CO2 flux parameters. TRACE uses an online full‐physics mesoscale atmospheric model and assimilates observations serially in a coupled atmosphere–carbon ensemble Kalman filter. The data assimilation system was tested in a series of observing system simulation experiments using in situ observations for a regional domain over North America in summer. Under ideal conditions with known prior flux parameter error covariances, TRACE reduced the error in domain‐integrated monthly CO2 fluxes by about 97% relative to the prior flux errors. In a more realistic scenario with unknown prior flux error statistics, the corresponding relative error reductions ranged from 80.6% to 88.5% depending on the specification of prior flux parameter error correlations. For regionally integrated fluxes on a spatial scale of 106 km2, the sum of absolute errors was reduced by 34.5%–50.9% relative to the prior flux errors. Moreover, TRACE produced posterior uncertainty estimates that were consistent with the true errors. These initial experiments show that the ESSPE approach in TRACE provides a promising method for advancing CO2 inversion techniques.

the atmosphere (also known as CO 2 fluxes) are relatively well understood on a global scale (Friedlingstein et al., 2021), there are major knowledge gaps at regional scales, which limit our ability to predict future climate change (e.g., Friedlingstein et al., 2014). Terrestrial ecosystems have removed about one-third of human-produced atmospheric CO 2 , but whether this terrestrial carbon sink is mainly driven by longer growing seasons (Piao et al., 2007), the CO 2 fertilization effect (C. Chen et al., 2022;Schimel et al., 2015), increased soil nutrient availability due to elevated nitrogen deposition (Lu et al., 2021;Reay et al., 2008), regrowth of previously logged forests in mid-to-high latitudes (Pugh et al., 2019), or other factors remains an open question, in part because the spatial distribution of this sink is not known well. To address this and other important open issues in carbon cycle science, accurate estimates of regional-scale CO 2 fluxes are needed.
It is generally not possible to directly measure surface CO 2 fluxes at subcontinental to continental scales partially due to their strong heterogeneity. Observations using the eddy covariance technique are typically representative of fluxes within a few hundred meters of the measurements (Chu et al., 2021), and upscaling of such observations can lead to significant biases (Ran et al., 2016). Correctly modeling the net flux of CO 2 between the terrestrial biosphere and the atmosphere, known as net ecosystem exchange (NEE), is also challenging because NEE is the residual of two typically large fluxes, uptake through gross primary production and release through ecosystem respiration. A relatively small error in simulating either of these fluxes can lead to significant systematic errors in NEE on monthly or longer time scales. Seiler et al. (2022) found that biosphere models do a reasonably good job in simulating the terrestrial carbon sink, but nevertheless leave room for improvements in terms of biases and regional-scale processes.
An alternative approach for deriving surface CO 2 fluxes is to use atmospheric observations of CO 2 concentrations. The atmosphere is a natural integrator of CO 2 fluxes through atmospheric mixing and transport, and variations in CO 2 concentrations higher up in the atmosphere therefore tend to be representative of CO 2 fluxes over larger regions. In this approach, known as inverse modeling, surface CO 2 fluxes are inferred from variations in atmospheric CO 2 using an atmospheric transport model. Early inversion studies focused primarily on the global-scale carbon cycle using CO 2 concentration measurements sampled in locations remote from anthropogenic and natural sources and sinks (e.g., Tans et al., 1990). Subsequent studies found a wide spread in inversion results on continental to subcontinental scales (Gurney et al., 2002;Peylin et al., 2013), which was partly attributed to differences in the modeling of atmospheric CO 2 transport. The expansion of the in situ CO 2 network and availability of new remote sensing observations of atmospheric CO 2 have motivated the development of regional inversion systems (e.g., Gourdji et al., 2010;Lauvaux et al., 2012;Monteil & Scholze, 2021;Schuh et al., 2010;Wesloh et al., 2020). Compared with global systems, regional inversions can represent the CO 2 fluxes and atmospheric transport at higher resolutions at feasible computational costs. However, even in regions with relatively well-developed observational networks such as Europe, large discrepancies in subcontinental CO 2 fluxes still remain in regional inversions (Monteil et al., 2020).
In light of the Paris Agreement, there has recently been a growing interest in using inversion methods to support national greenhouse gas monitoring and reporting (IPCC, 2019;Janssens-Maenhout et al., 2020). The number of CO 2 observations is expected to increase severalfold in the next decade to address this need. For example, a constellation of operational satellites code-named the Copernicus Anthropogenic Carbon Dioxide Monitoring mission (CO2M for short) will deliver data on column-averaged dry air mole fraction of CO 2 (XCO 2 ) with a swath width of about 250 km and at 2 × 2 km spatial resolution (Meijer, 2019). Using such satellite observations is computationally prohibitive for most common inversion schemes, especially for inversions that use Lagrangian backward models to compute the sensitivity of observations to surface fluxes, as the satellite retrievals are not only numerous but also sensitive to the whole atmospheric column. The computational costs can be reduced using techniques that for example, calculate the observational sensitivity for only a subset of the soundings (e.g., Roten et al., 2021), but such methods are usually associated with their own limitations, such as not being applicable to broader regions or longer time periods.
Another challenge associated with XCO 2 retrievals from satellites is how to account for correlated random and systematic errors in the spatially dense soundings. Correlated errors in XCO 2 can arise due to errors in both the instruments and retrieval process (e.g., Baker et al., 2022), as well as due to errors in the simulated atmospheric transport in the inversion. In most inversion schemes, atmospheric transport errors would be prescribed in the observation error covariance matrix, which is typically assumed to be diagonal (i.e., uncorrelated observation errors) to simplify the calculations. Chevallier et al. (2010) showed that neglecting correlated atmospheric 10.1029/2022MS003208 3 of 20 transport can cause significant biases in the inverted fluxes when assimilating XCO 2 observations, even when using techniques such as observation thinning.
The significantly improved availability of CO 2 observations in the near future and new applications of CO 2 inversions underscore the need to develop inversion schemes that are capable of (a) assimilating a large number of observations, including column-averaged observations from satellites; and (b) accounting for correlations in the observation and model errors, for instance correlated atmospheric transport errors. One promising approach is to couple a CO 2 inversion scheme with an atmospheric data assimilation system and jointly optimize for surface CO 2 fluxes (or parameters controlling the fluxes) and atmospheric variables including atmospheric CO 2 concentrations. In this dual state and parameter estimation formulation, CO 2 concentrations are considered as part of the state, and CO 2 fluxes are treated as parameters. Uncertainties in modeled CO 2 concentrations, which arise partly due to atmospheric transport errors, would then be included in the prior error covariance matrix, which, unlike the observation error covariance matrix, is allowed to be non-diagonal. Kang et al. (2012) implemented a dual-state CO 2 flux inversion system using the local ensemble transform Kalman filter and a simplified global atmospheric transport model, and found that the system could accurately estimate the synthetic true surface CO 2 fluxes using idealized satellite XCO 2 retrievals. Including the error covariances between winds and atmospheric CO 2 concentrations was found to be beneficial, while error covariances between CO 2 concentrations or surface fluxes and other meteorological variables degraded the results (Kang et al., 2011). Other examples of global dual-state CO 2 inversion systems include the Tan-Tracker system (Tian et al., 2014) and a system based on the local ensemble transform Kalman filter (Y. Liu et al., 2019). At regional scale, Peng et al. (2015) and Kou et al. (2017) tested two systems based on the offline Community Multiscale Air Quality (CMAQ) transport model, both at a horizontal resolution of 64 km × 64 km. There is a need to further test and develop dual-state CO 2 inversion algorithms, especially using online atmospheric transport models, which would allow for higher-resolution simulations of the atmospheric CO 2 transport. This paper presents the TRACE Regional Atmosphere-Carbon Ensemble (TRACE) system, a regional coupled atmosphere-carbon data assimilation system for CO 2 flux inversion. TRACE uses the ensemble-based simultaneous state and parameter estimation (ESSPE) approach to jointly optimize atmospheric CO 2 mole fractions and parameters controlling surface CO 2 fluxes. While the inclusion of atmospheric CO 2 mole fractions in the state vector makes it possible to reduce the uncertainty in the initial atmospheric CO 2 field (e.g., Tian et al., 2014), the main motivation for the dual state and parameter estimation approach here is to make it computationally feasible to assimilate a large number of observations from for example, CO 2 -observing satellites. TRACE uses an online mesoscale numerical weather prediction model to dynamically downscale meteorological variables from a global data set, and additionally includes an atmospheric data assimilation component, which makes it possible to simulate and reduce atmospheric transport errors by assimilating meteorological observations (see H. W. Chen et al., 2019). The objective of this study is to introduce the data assimilation methods used by TRACE and comprehensively test the system in a series of controlled experiments with synthetic in situ observations. More advanced features, such as assimilation of XCO 2 observations from satellites and inclusion of atmospheric transport errors, are left for future studies.
Section 2 provides a general overview of the TRACE framework, followed by detailed descriptions of the major components in the subsections. Section 3 explains the specific experiments performed in this study, and Section 4 presents the results. Finally, Section 5 summarizes and concludes the paper. Figure 1 shows a schematic of the key components of TRACE. Central to the modeling framework is the regional atmospheric model (Section 2.2), which simulates the mixing and transport of CO 2 in the atmosphere. The atmospheric model downscales meteorological fields from a global data set (Section 3.1) to the mesoscale (spatial scales of ∼1-100 km). Atmospheric CO 2 mole fractions are provided to the model as initial and lateral boundary conditions, and surface CO 2 fluxes as lower boundary conditions.

The TRACE Framework
In the inversion procedure, a set of parameters controlling the CO 2 surface fluxes are optimized to reduce the mismatches between simulated CO 2 mole fractions and atmospheric CO 2 observations. The parameters are currently simple scaling factors on prior CO 2 fluxes (Section 2.1), but in principle it is also possible to optimize model parameters in for example, a vegetation model, similar to a carbon cycle data assimilation system approach (Scholze et al., 2007). The parameter values are estimated at the pixel scale on the same spatial resolution as the atmospheric model.
The inversion procedure is divided into two steps: forward modeling and assimilation. The purpose of the forward modeling step is to derive the sensitivities of atmospheric observations to CO 2 flux parameters, which is achieved by running an ensemble of forward simulations with perturbed flux parameters using the atmospheric model. The flux parameter perturbations are randomly drawn according to the prescribed prior flux parameter error covariances. In the TRACE framework, it is moreover possible to include perturbed meteorological fields in the ensemble members to simulate and constrain atmospheric transport errors (H. W. Chen et al., 2019). This feature of the coupled atmosphere-carbon data assimilation system will be investigated in future studies.
The data assimilation component is based on a deterministic square-root ensemble Kalman filter (EnKF), see Section 2.3. During the assimilation step, the perturbed flux parameters and corresponding simulated atmospheric CO 2 mole fractions from the ensemble members are fed to the data assimilation component as prior estimates. The simulated atmospheric CO 2 mole fractions are translated to observation equivalents using an observation operator, which for in situ observations corresponds to straightforward linear interpolations in time and space, and are then compared with the observations to produce a vector of data-model differences, also known as the innovation. The data assimilation algorithm computes a correction to the prior CO 2 flux parameters and mole fractions-the analysis increment-by applying a weighting matrix to the innovation, taking into account the uncertainties and error covariances of the prior estimates and observations. In the ESSPE approach, the CO 2 Figure 1. Simplified schematic of the TRACE system. The data flow shows how the prior CO 2 fluxes in the top left are propagated through the modeling chain and combined with other data to derive the posterior fluxes in the bottom right. Parallelograms denote data (gray parallelograms were not included in this study) and orange rectangles denote model components. The dashed lines indicate how the ensemble of posterior atmospheric fields, including atmospheric CO 2 mole fractions, and posterior flux parameters are passed to the next assimilation cycle as the new prior. Some modeling steps are not shown for simplicity, such as covariance localization and inflation, reinitialization of meteorological fields, and forecasting of flux parameters. mole fractions are jointly updated with the CO 2 flux parameters to reflect the updated parameter values, and thus there is no need to re-run the atmospheric model with the new parameter values. In addition to updating the ensemble mean estimates to be closer to observations, the ensemble deviations are also adjusted to reflect the reduced uncertainty after assimilating observations. The optimized CO 2 flux parameters and mole fractions are then fed back to the atmospheric model, and the forward modeling and assimilation steps are repeated for the next data assimilation cycle.

Flux Model
The flux model simulates the CO 2 fluxes based on a set of input data and flux parameters. Currently, the space-and time-varying total CO 2 fluxes F(x, y, t) are calculated based on prior fluxes from different data sets, which are scaled by flux parameters for the different flux components: where s bio , s ocn , s ff , and s fire are scaling factor parameters for the prior terrestrial biogenic fluxes (F bio ), oceanic fluxes (F ocn ), fossil fuel emissions (F ff ), and wildfire emissions (F fire ), respectively. Section 3.1 provides more information about the specific prior fluxes used in this study.
In this paper, we considered only s bio and s ocn for optimization. To simplify the interpretation of the results we did not include errors in the fossil fuel or wildfire emissions; these could be included in future studies in a straightforward manner. Three different error correlation functions were tested for the prior flux parameter errors: (a) ecoregion-specific parameters, with fully correlated parameter errors within the same ecoregions, and no correlations between ecoregions; (b) distance-based error correlations that decrease with distance; and (c) a hybrid of the previous two approaches, with distance-based error correlations within ecoregions, and an additional decay factor that decreases the error correlations between ecoregions. The ecoregions in the domain presented in this study are shown in Figure 2. Figure 3a shows how the error correlations vary with distance and between ecoregions for the three correlation functions.
The first ecoregion-based correlation function is based on the consideration that many vegetation models calculate CO 2 fluxes based on parameters that are assumed to be specific to certain ecoregions or plant functional types. The ecoregions in this study are based on the ecosystem classification by Olson et al. (1985) and were further simplified to reduce the number of ecoregions (Text S1 in Supporting Information S1). The ocean basins were treated in a similar way to ecoregions. We ignored CO 2 flux parameters for the Desert and Tundra ecoregions and inland water bodies following Jacobson et al. (2020). An example of the ecoregion-based flux error correlations for a grid point in the Cool Mixed ecoregion is shown in Figure 3b.
The ecoregion-based correlation function makes use of prior information about the distribution of ecoregions, but may overestimate the error correlations between spatially distant grid points because they happen to share the same ecoregion. Moreover, this approach creates sharp contrasts in the spatial correlations (see Figure 3b), which could add unphysical updates to the posterior fields and make the inversion more sensitive to atmospheric transport errors.
The second distance-based correlation function considers only the distance between grid points, as shown in Figure 3c. The error correlations r decay with distance d according to a Gaussian function: where L is the correlation length. We also tested an exponentially decaying correlation function but found that it led to less smooth updates to the posterior fields due to the fast correlation decay at short distances (see Figure S1 in Supporting Information S1), and hence only show results for the Gaussian function. In this study, the correlation length was set to 500 km for land ecoregions and 3,000 km for ocean basins. The correlation between ecoregion and ocean flux parameters was assumed to be zero. Correlations below 0.01 were set to zero to reduce long-distance correlations. While the distance-based correlation function results in smoother updates to the posterior fields, the completely isotropic correlations could potentially result in unphysical error correlations between different ecoregions.
The hybrid correlation function is an attempt to overcome deficiencies in the previous two approaches. It is similar to the distance-based correlation function, except for the addition of a decay factor γ when calculating the error correlations between grid points belonging to different ecoregions. In the hybrid formulation, the correlation between two grid points of different ecoregions is calculated as: where L d is the decay length, and d eco is the average of two distances, d 1 and d 2 . For two grid points x 1 and x 2 belonging to two different ecoregions eco 1 and eco 2 , d 1 is the shortest distance between x 1 and a grid point belonging to eco 2 , and similarly for d 2 but for x 2 and eco 1 . Here the hybrid correlation function was tested with L = 600 km and L d = 300 km. The correlation length L was set to a higher value than for the distance-based correlation function on the basis that we may expect longer-range correlations for grid points that share the same ecoregion. Figure 3d shows that the hybrid-based correlation function retains some information about the ecoregions, while simultaneously creating smoother correlation gradients between ecoregions and reducing far-distant correlations. The hybrid approach is flexible in the sense that it can reproduce both the ecoregion-and distance-based correlation functions: letting L d approach infinity reduces Equations 2 to 3, while letting L approach infinity and L d approach zero yields the ecoregion-based correlation function.
In this study, the flux parameters were assumed to have a prior uncertainty (one standard deviation) of 0.8 (unitless) for land ecosystems (s bio ) and 0.4 (unitless) for ocean flux parameters (s ocn ) following CarbonTracker CT2019B (Jacobson et al., 2020). The CO 2 flux parameters were treated as constants and were forecast using a simple persistence model. The ensemble mean parameter values were not relaxed back to the prior or initial parameter values, as is done in for example, CT2019B. Such relaxation techniques could be useful to avoid longterm biases in parameters that are not well constrained by observations, and will be considered and expanded upon in future studies.

Atmospheric Model
The atmospheric model in the TRACE system is the Weather Research and Forecasting (WRF) model coupled with Chemistry (WRF-Chem; Grell et al., 2005) version 3.6.1. WRF-Chem is a regional chemical transport model based on the mesoscale WRF model (Skamarock et al., 2008). The regional atmospheric model is used for primarily three purposes: (a) to compute the sensitivities of atmospheric CO 2 observations to flux parameters; (b) to dynamically downscale meteorology from a global data set to the mesoscale, which can improve the modeling of atmospheric CO 2 transport (Agustí-Panareda et al., 2019); and (c) to act as the dynamic model for atmospheric variables and CO 2 mole fractions, making it possible to simulate the error growth in these variables due to uncertainties in meteorological fields and model physics. The WRF-Chem version used in TRACE has been modified to allow for multiple sets of CO 2 fluxes and tracers within one single run, making it computationally efficient to perform ensemble runs with perturbed CO 2 flux parameters. Atmospheric CO 2 mole fractions are provided from a coarse-resolution global model as lateral boundary conditions for the limited-area model using a mass-conserving interpolation algorithm (Butler et al., 2020).
In this study, WRF-Chem was set up over North America (see Figure 2) on a spatial resolution of 27 km and with 60 vertical levels. The time step was set to 60 s. Atmospheric CO 2 was modeled as an inert passive tracer with no chemical reactions and no feedback to meteorology. The model configuration and physical parameterization schemes followed the setup in H. W. Chen et al. (2019), including a positive definite sixth-order diffusion scheme, the Noah land-surface model (F. Chen & Dudhia, 2001), the Mellor-Yamada Nakanishi and Niino Level 2.5 boundary layer scheme (Nakanishi & Niino, 2006), and the Kain-Fritsch convective scheme (Kain, 2004) with no parameterized convective transport for tracers.

Data Assimilation
The data assimilation system is based on the EnKF algorithm implemented in the Pennsylvania State University (PSU) WRF EnKF system (Zhang et al., 2006). This section gives a brief overview of the general data assimilation algorithm, and Section 2.4 presents the specifics of the dual state and flux parameter estimation approach for CO 2 inversions.
The Kalman filter is formulated in a Bayesian framework, where the optimal estimate of the state of a dynamical system is derived by combining information from a dynamic model and observations, taking their respective uncertainties into account. First, a prediction is made using the dynamic model. The prediction step is followed by the analysis step, during which the state vector x containing state variables from the dynamic model is updated in the Kalman filter according to where x a is the updated posterior (or analysis) state vector after combining the prior model estimate and observations, x b is the prior (or background) state vector from the model prediction, y o is a vector containing all observations, H is the linearized observation operator relating model variables to observations, and K is the Kalman gain matrix that minimizes the expected variance. For a linear system with Gaussian errors, the Kalman gain takes the form where P b and R are the error covariance matrices for the prior (background) and observations, respectively. The Kalman gain determines the weights of the information provided by the prior and observations, and spreads the influence of observations in model space according to the prior error covariances. At the end of the analysis step, the prior error covariance matrix is updated to reflect the reduced uncertainty from assimilating the observations according to: where P a is the posterior (analysis) error covariance matrix, and I is the identity matrix.
The EnKF is a Monte Carlo approximation of the Kalman filter that is computationally feasible for high-dimensional systems with very large state vectors. In the EnKF, the prediction step is carried out through an ensemble of predictions using the full (potentially nonlinear) dynamic model. The prior error covariance matrix is replaced by the ensemble approximation where N is the number of ensemble members and X′ is a matrix whose columns are state vectors of the ensemble perturbations (deviations from the ensemble mean).
It is often computationally challenging to invert the matrix in parentheses in Equation 5, which has the size p × p, where p is the number of observations. The data assimilation component in TRACE uses a deterministic square-root formulation of the EnKF to make it possible to assimilate a large number of observations. In this formulation, observations are assimilated serially under the assumption that R is diagonal, thus reducing the matrix inversion in Equation 5 to a scalar division. The ensemble mean state ̄ is updated iteratively for each observation using Equation 4, replacing all x in the equation with ̄ . In each iteration, the ensemble perturbations are also updated to correspond to the new posterior error covariances after assimilating the observations: where x′ a and x′ b are the posterior and prior ensemble perturbations, respectively, and ̃ is the reduced Kalman gain:̃= Because observations are assimilated one at a time, α, R, and HP b H ⊤ are scalars in Equation 9. It is also worth noting that the serial EnKF allows the use of the full observation operator  rather than the linearized version H. For our purposes the observation operator is already linear (see Section 2.4 for more information), thus it makes no difference whether we use  or H. After the analysis step, the ensemble members are propagated to the next analysis time using the dynamic model, and the whole cycle repeats.
In EnKF applications, the dimension of the state vector is typically much larger than the ensemble size, which can lead to rank deficiency and spurious correlations between unrelated state vector elements. This issue is often addressed by localizing the impact of observations to remove artificial correlations far away from the observations. Another potential issue is that the posterior errors are underestimated due to for example, neglected error sources, which can eventually lead to a collapse of the ensemble spread, also known as filter divergence. In this condition the EnKF will ignore the observations in Equation 4 because P b , and therefore K, approaches a zero matrix. To avoid filter divergence, the ensemble spread is often inflated before or after the analysis step. The specific localization and inflation methods used in this study are detailed in Section 2.5. For more information about the EnKF, we refer to Houtekamer and Zhang (2016).

Ensemble-Based Simultaneous State and Parameter Estimation
The EnKF algorithm in TRACE has been well tested for numerical weather prediction (e.g., Zhang et al., 2006;Zhang et al., 2019). To extend the data assimilation algorithm for CO 2 flux inversion purposes, an ESSPE approach was used, whereby the state vector was expanded to include both atmospheric CO 2 mole fractions at every grid point and model level, as well as parameters controlling surface CO 2 fluxes (see Section 2.1). The act of including parameters in the state vector is called state augmentation (e.g., Annan, 2005). Here, the underlying assumption is that there are persistent errors in the flux parameters that vary slowly in time. The flux parameters differ from state variables from the dynamic model in two important aspects: (a) there is no prognostic equation for the parameters in the dynamic model, and (b) the parameters are not physical quantities that can be directly observed. In conventional inversions that include only fluxes or flux parameters in their state vectors, the fluxes or parameters are related to observations through an atmospheric transport model as part of the observation operator H. For the ESSPE approach, however, the links between flux parameters and atmospheric CO 2 mole fractions are instead expressed in the prior error covariance matrix P b , which is estimated from the ensemble members (see Equation 7). Consequently, the observation operator relating modeled CO 2 mole fractions to in situ CO 2 observations becomes a simple interpolation operation in time and space to the observation points.
An additional advantage of the ESSPE approach is that the prior error covariance matrix P b can contain other error sources for simulated atmospheric CO 2 mole fractions, including atmospheric transport errors and uncertainties in CO 2 lateral boundary conditions in regional inversions. Such errors can exhibit correlations over long distances and across time (e.g., H. W. Chen et al., 2019;Feng et al., 2019;McNorton et al., 2020). However, conventional inversions typically prescribe these errors in the observation error matrix R, which is often assumed to be diagonal (uncorrelated errors) to make the calculations tractable. In the ESSPE approach such errors would instead be included by adding additional perturbations to the ensemble members, and the off-diagonal covariances would be preserved in P b .
The classical EnKF operates on the state vector one assimilation time step at a time. There are variations of the EnKF called ensemble Kalman smoothers that consider the state and/or observations over an assimilation window. In this study we performed CO 2 flux inversions in TRACE using both a filter with an assimilation time step of 3 hr and no assimilation window, as well as a fixed-lag smoother with a window size of 1 week. For the smoother, TRACE assimilated 3-hourly observations over 1-week windows and optimized 3-hourly CO 2 mole fractions and weekly flux parameters over the same time windows as the observations, with no overlaps between the windows.

Covariance Localization and Inflation
A horizontal localization scheme was applied to in situ CO 2 observations with the shape of the Gaspari and Cohn (1999) fifth-order polynomial, which has a similar shape as a Gaussian function but with a faster convergence to zero at the tails. The radius of influence was set to 5,400 km; the relatively large radius of influence was found to be beneficial to optimize distant flux parameters. The same localization function was applied when updating both CO 2 mole fractions and flux parameters, and no vertical localization was used for CO 2 mole fractions.
When optimizing ecoregion-specific flux parameters, the distance-based covariance localization could lead to inconsistent updates to parameter values within the same ecoregion. Thus, when using the ecoregion-based correlation function, the parameter values within each ecoregion were averaged after each data assimilation procedure to obtain a single parameter value for every ecoregion.
To avoid filter divergence, the ensemble spread (including both CO 2 mole fractions and flux parameters) was inflated after each analysis step by relaxing the posterior ensemble perturbations to their prior values following Zhang et al. (2004): where ′ relaxed is the new posterior ensemble perturbation for a particular ensemble member and α is the relaxation coefficient. In this study α was set to 0.2.
Even with this inflation scheme, the estimated errors in the parameters would inevitably decrease after each assimilation cycle as there is no error growth in the parameters due to the lack of error growth in the dynamic model for the flux parameters. Hence, the parameters were also partially relaxed to their prescribed initial perturbations after each analysis step (except for the last analysis step) in preparation for the next cycle according to ′ relaxed = sgn where sgn(x) is the element-wise sign function (−1 if x is negative, 0 if x is 0, and 1 if x is positive), abs(x) is the element-wise absolute value function, and β is the relaxation coefficient for parameter values. The value of β was set to 0.1 in our experiments. s′ is a vector containing perturbed parameter values, with ′ relaxed , ′ , and s′ i denoting the new relaxed parameter perturbations, the posterior parameter perturbations, and the initial parameter perturbations prescribed at the start of the experiments, respectively. Thus, at a minimum, the parameters will have an uncertainty of 10% of the initial uncertainty. The parameter inflation plays a critical role in preserving the links between the flux parameters and CO 2 mole fractions in the off-diagonals of the prior error covariance matrix estimated from the ensemble members (Equation 7).
In Equation 11, the sgn(x) and abs(x) functions are included to account for cases with nonlinearities, for example, due to ensemble members having different atmospheric transport, which can result in corresponding elements in s′ a and s′ i having opposite signs. If that happens, the relaxation scheme could cause the ensemble spread to decrease instead of increase. In this study we considered only cases with perfect transport, in which case Equation 11 takes the same form as Equation 10.

Observing System Simulation Experiments
A series of observing system simulation experiments (OSSEs) were performed to thoroughly test the ESSPE approach in TRACE for regional CO 2 inversions. The experiments focused on constraining biogenic and oceanic CO 2 fluxes in North America during the month of July in 2016. We chose this case because first, the biogenic activity is strongest in summer, and second, it coincides with the first Atmospheric Carbon and Transport-America (ACT-America) airborne field campaign (Davis et al., 2021), which can provide additional observations for assimilation and verification in future studies using real data.
The OSSEs were performed by first running a reference simulation with one set of CO 2 fluxes, which was taken as the truth. Synthetic observations were then sampled with noise from the truth simulation according to a realistic network of ground-based CO 2 monitoring sites (see Figure 4). Next, we assimilated the synthetic observations in data assimilation experiments to assess how well the TRACE system could reconstruct the true fluxes. These experiments used different prior CO 2 fluxes than the truth simulation but were otherwise identical to the reference. Further information about the data and experimental setup is provided in the following subsections.

Atmospheric and Flux Data
The regional atmospheric model was driven by meteorological initial and boundary conditions from the ERA-Interim reanalysis (Dee et al., 2011). The ERA-Interim data had a spatial resolution of about 81 km, 60 vertical levels, and 6-hourly time resolution. Initial and lateral boundary conditions in CO 2 mole fractions were taken from the global CarbonTracker CT2019B data set (Jacobson et al., 2020), which has a spatial resolution of 3° × 2° longitude-latitude and 25 vertical levels. The initial conditions in atmospheric CO 2 mole frac tions were created by first initializing WRF-Chem with mole fractions from CT2019B (linearly interpolated to the WRF-Chem grid), and then running the model for two months to allow the mesoscale model to simulate high-resolution gradients in atmospheric CO 2 . During the two-month spin-up, the meteorological fields were reinitialized from ERA-Interim every five days with an initial 12 hr spin-up each time.
The prior CO 2 fluxes came from the prior fluxes in CT2019B with a spatial resolution of 1° × 1°. The biogenic fluxes were originally from a Carnegie-Ames Stanford Approach (CASA; Potter et al., 1993) biogeochemical model simulation called CASA GFED 4.1s, and the oceanic fluxes were based on the Takahashi et al. (2009) climatology. These monthly-mean fluxes had been postprocessed by the CarbonTracker team and downscaled to 3-hourly time resolution, as documented on the CT2019B website (http://dx.doi.org/10.25925/20201008). We spatially interpolated the CT2019B prior fluxes to the WRF-Chem grid using a first-order conservative regridding scheme. Because we chose to exclude the uncertainties in fossil fuel and wildfire emissions, these flux components would have no effect in our perfect transport OSSEs, and were thus set to zero. In the final set of experiments, we used different fluxes for the truth simulation based on an alternative CASA simulation (CASA GFED_CMS) for biogenic fluxes and Ocean Inversion Fluxes (Jacobson et al., 2007) for oceanic fluxes, which are also used as prior fluxes in CT2019B.

Observations
TRACE was tested with synthetic in situ CO 2 observations, the locations of which are shown in Figure 4. We focused primarily on NOAA and Environment Canada towers (see Table S1 in Supporting Information S1 for a full list) that were available during the study period. The observations were sampled from the truth simulation as 3-hourly instantaneous CO 2 mole fractions at the locations and measurement heights of the real measurements. To simulate observation errors due to representativeness errors, errors in the observation operator, and measurement errors, we added random errors to the observation values drawn from a Gaussian distribution with a mean of 0 and standard deviation of 1 ppm.
Many inversion systems use predefined times for when observations are considered for assimilation, for example, using only afternoon values for observations within the planetary boundary layer (PBL) to ensure well-mixed conditions. Here we used an alternative approach where instead of specifying the times a priori, we considered all observations and filtered out those that were deemed to not be representative of fluxes at the scales of interest. Specifically, we kept the observations that met at least one of the following criteria in summer: (a) the modeled PBL height at the observation site was 1,000 m or higher; (b) the observation was within the marine boundary layer; or (c) the observation was at least 100 m above the PBL in the model. The 100 m margin was included to account for possible errors in modeled PBL heights in real cases.
An observation error of 4 ppm was assigned to the in situ CO 2 observations in the observation error covariance matrix when assimilating the observations. The observation error was set to a larger value than the true observation errors of 1 ppm to take into account that the specification of prior flux parameter error covariances is imperfect, which can be considered as part of representativeness errors in the model. Erroneous prior flux parameter error covariances can lead to incorrect upscaling of the flux corrections, and a larger observation error gives the inversions more room to avoid overfitting to observations in such cases.

Experiments
A total of five OSSEs were performed, which are listed in Table 1. In the first two experiments (ECO and DIST), the true fluxes were generated by perturbing the prior fluxes with random scaling factors (80% uncertainty on land parameters and 40% on ocean parameters, see Section 2.1) sampled from the same prior flux parameter error covariance matrix as the ensemble members. In other words, we assumed that the true flux error covariances were known in these experiments. This assumption is highly over-optimistic; the purpose of these experiments was to test the data assimilation system under ideal conditions, which is useful to understand the "intrinsic" constraint on CO 2 flux parameters given nearly perfect conditions (analogous to intrinsic predictability, see e.g., Melhauser & Zhang, 2012).
For the last three ALT experiments, we used alternative datasets for the biogenic and oceanic fluxes (see Section 3.1) in the truth simulation, reflecting a more realistic scenario where the true flux error statistics are unknown. The inversions used the same prior fluxes as before. Three inversions with different prior flux param-

Experiment
True The data assimilation experiments were started on 3 June 2016 at 0000 UTC to allow the system to spin up for 4 weeks. First, we ran the atmospheric model for 2 weeks without assimilating observations to establish the initial sensitivities of atmospheric CO 2 mole fractions to the flux parameters. Next, in situ CO 2 observations were assimilated to jointly constrain the atmospheric CO 2 mole fractions and CO 2 flux parameters. The first 2 weeks of assimilation were discarded to account for spin-up. Because most experiments used an assimilation window of 1 week, the experiments ended on 5 August 2016, 9 full weeks after the start date. Posterior fluxes beyond July 2016 were discarded in the analyses to focus on the month of July.
We evaluated the performances of the inversions based on the errors in spatially integrated posterior fluxes. The error is defined as the difference between the posterior or prior fluxes and the true fluxes from the reference simulation. We also examined the estimated errors from the ensemble spread, which we will refer to as uncertainties to differentiate them from the true errors. For the errors (and uncertainties) we calculated the relative error (uncertainty) reduction, defined as 1 minus the ratio of the posterior and prior absolute errors (uncertainties). The closer the relative error reduction is to 100%, the better the inversion performed, while a negative error reduction indicates that the inversion degraded the prior fluxes.

Idealized Inversions With Known Prior Flux Parameter Error Covariances
In the first ECO experiment, the true NEE fluxes were generated by perturbing the prior NEE with ecoregion-based scaling factors (see Section 2.1 for details). Figure 5 shows the true scaling factors for six selected ecoregions along with the estimated scaling factors from the inversion. The prior parameter values were initially 1 for all ecoregions, and this was also the value for the estimated parameters during the first 2 weeks when no observations were assimilated. After 2 weeks the system started to assimilate in situ CO 2 observations to simultaneously optimize the atmospheric CO 2 state and CO 2 flux parameters, which brought the estimated scaling factors closer to the true values. For some ecoregions, for example, Warm Crops and Cool Crops (Figures 5a and 5b), the inversion initially slightly overcorrected the flux parameters, and required one additional assimilation cycle for adjustments. This overcorrection is more likely to occur at the beginning of the data assimilation window due to the large prior uncertainties and model-data mismatches before any observations have been assimilated. The first 2 weeks of assimilation were therefore discarded as spin-up.
In addition to improving the ensemble mean estimates of the flux parameters, the posterior uncertainties in the parameters were also continually updated after each assimilation cycle, as shown by the shadings in Figure 5.
In ecoregions with flux parameters that were relatively well constrained by the atmospheric observations, for example, Warm Crops (Figure 5a), Cool Crops (Figure 5b), and Warm Forest/Field (Figure 5d), the posterior uncertainties were also strongly reduced.
The flux parameter for the Tropical ecoregion in Figure 5f converged slower to the truth than the other flux parameters due to the low observational sensitivity to fluxes in the remote Tropical ecoregion. In the ESSPE approach, it is possible to constrain distant fluxes even with a relatively short 1-week assimilation window because (a) the errors in the flux parameters are assumed to be strongly correlated in time (here the parameters were treated as constants); and (b) the prior error covariances, which contain the link between the CO 2 mole fractions and fluxes, are cycled between assimilation windows. Thus, later observations can constrain past parameter values as air from distant regions arrives at the observation sites. This aspect is shown clearly by the experiment using a filter without an assimilation window (orange line in Figure 5), which was still able to estimate flux parameters in areas remote from observations such as the Tropical ecoregion. The filter inversion exhibited stronger fluctuations in estimated flux parameters at the start of the experiment, which were reduced when using a 1-week assimilation window because some of the random errors were averaged out over the assimilation window. We therefore show results from inversions that used the assimilation window setup in the subsequent analyses. Figures 6a-6c show the spatial distribution of the monthly NEE from the ECO experiment for the prior fluxes, true fluxes, and posterior fluxes after assimilating observations, respectively. The fluxes were integrated over the month of July 2016. When both the true fluxes and the perturbed fluxes in the ensemble members were derived using ecoregion-based scaling factors, the system was able to accurately retrieve the true fluxes (compare Figures 6c and 6b). The difference in monthly domain-integrated NEE between the posterior and truth is −3.2 Tg C month −1 , compared with −134.1 Tg C month −1 for the prior and truth, corresponding to a 96.9% relative error reduction. Although the ensemble-mean scaling factor for the Tropical ecoregion did not reach the true value (Figure 5f), the NEE in this ecoregion was relatively small and therefore had a negligible effect on the error in domain-integrated NEE. In terms of spatial correlation, the correlation between the posterior fluxes and the truth is r = 1.0, compared with r = 0.77 for the prior. The inversion using a filter setup showed similar results (not shown).
For the ECO experiment there were 15 flux parameters to optimize, one for each ecoregion or ocean basin. To test a more complex case, the flux parameters were instead assumed to vary per grid point with a distance-based error correlation (500 km correlation length for land fluxes) in the DIST experiment (see Section 2.4 for details). Figure 6e shows the true NEE for the DIST experiment and Figure 6f the posterior NEE after assimilating observations. There are some differences between the posterior and true fluxes, most noticeable in western Wyoming and Colorado where the inversion overestimated the strength of the CO 2 sink, likely due to the scarcity of observations around this region. Nevertheless, the posterior fluxes are considerably closer to the truth than the prior fluxes in terms of both spatial pattern (r = 0.87, compared with r = 0.68 for the prior fluxes) and domain-integrated NEE; the difference from the truth is −2.9 Tg C month −1 , compared with −134.1 Tg C month −1 for the prior fluxes, corresponding to a 97.9% relative error reduction. Repeating the inversion with an extensive theoretical in situ network-tower observations at 100 m above ground level spaced 540 km apart-reduced the erroneous differences in the posterior fluxes and increased the spatial correlation to r = 0.93 (not shown).
In addition to improving the CO 2 flux estimates, the ESSPE algorithm simultaneously reduced the errors in simulated atmospheric CO 2 mole fractions. In the ECO experiment, the domain-integrated root-mean-square error in XCO 2 (excluding 20 grid points around the edges of the domain) was reduced, on average, from 0.65 ppm in the simulation with prior fluxes, to 0.03 ppm after assimilating the CO 2 observations ( Figure S2a in Supporting Information S1). For the DIST experiment, the corresponding root-mean-square error in XCO 2 was reduced on average from 0.76 to 0.23 ppm ( Figure S2b in Supporting Information S1).

Inversions With Unknown Prior Flux Error Statistics
In the previous experiments, the true flux parameters and ensemble perturbation parameters were sampled from the same error covariance matrix, and with a sufficiently large ensemble the truth would lie in the subspace spanned by the ensemble members. In reality, the joint distribution of the true flux errors are unknown and largely uncertain. A new set of experiments were therefore carried out to test a more realistic scenario where the true flux errors and their statistics were not known a priori. The prior fluxes in these experiments were the same as in the previous experiments, while the true fluxes were obtained from different datasets than the prior fluxes.  that inversions using this type of prior flux parameter error covariances require a larger covariance inflation for the flux parameters. Overall, the posterior uncertainties from the ALT_DIST and ALT_HYBRID inversions correspond well to the true errors.
The time evolutions of daily domain-integrated NEEs are shown in Figure 8 for the prior, true, and posterior fluxes from the ALT experiments. All fluxes exhibit similar daily variations because they were downscaled to 3-hourly NEE using the same algorithm and driving data. The domain-integrated prior fluxes show an overall overestimated negative NEE throughout the period of interest. The negative NEE bias is substantially reduced in all three ALT inversions after assimilating the in situ CO 2 observations. Some negative biases persist in the posterior fluxes, particularly in the beginning and end of July. In the beginning of July, the biases in the posterior fluxes are mainly caused by small-scale errors in the Corn Belt region and negative NEE biases in the Southern United States and southeastern Mexico, while the biases at the end of July are mostly due to negative NEE biases in the Southern United States (not shown).

Regional Fluxes From Inversions
Next, we assessed the performance of the inversions in terms of constraining regional fluxes. A grid-scale comparison between the prior and posterior fluxes with the truth can be misleading because the prior and true fluxes have the same spatial resolution, while the posterior fluxes from the inversions introduced additional small-scale variations on the spatial resolution of the atmospheric model. Thus, metrics such as spatial correlation and root-mean-square error would tend to show unfavorable results for the high-resolution posterior fluxes. Instead, we divided the study domain into six regions of roughly the same area (5.3-5.9 × 10 6 km 2 ) and calculated the total NEE in each region. Figure 9a shows the regional NEEs for the prior, true, and posterior fluxes from the ALT experiments, and the corresponding domain-integrated NEEs are shown in Figure 9b. The six regions are displayed in Figure 9c together with the spatial distribution of the true fluxes and locations of observations. As seen from Figure 9a, the inversions were generally able to reduce the errors in the regional fluxes compared with the prior. Two exceptions are the NW and N regions, where the combination of positive and negative NEE biases in the prior result in a regionally integrated prior NEE that is close to the truth. In the N region, only the ALT_HYBRID inversion produced posterior fluxes with a smaller error (12.2 Tg C month −1 ) than the prior fluxes (−16.6 Tg C month −1 ). The overall performance in capturing the regional NEEs was quantified by the sum of absolute errors in regional fluxes: 157.7 Tg C month −1 for the prior fluxes, 103.3 Tg C month −1 for ALT_ECO, 92.6 Tg C month −1 for ALT_DIST, and 77.5 Tg C month −1 for ALT_HYBRID. Thus, in terms of regional fluxes, all inversions improved upon the prior fluxes, with the largest relative error reductions in the ALT_HYBRID inversion (50.9%), followed by ALT_DIST (41.3%) and ALT_ECO (34.5%).
For the domain-wide NEE, the posterior fluxes from the inversions were substantially closer to the true values than the prior fluxes in all experiments (Figure 9b). The relative error reductions in the domain-integrated NEE are 80.6% for the ALT_ECO experiment, 88.5% for ALT_DIST, and 84.4% for ALT_HYBRID. Although the ALT_HYBRID inversion has the smallest error in regional fluxes, the ALT_DIST inversion performed marginally better in terms of estimating the domain-wide NEE.
Figures 9a and 9b also show the estimated uncertainties (two standard deviations) in the prior and posterior fluxes as error bars. The prior uncertainties in regional NEE vary slightly between the three ALT experiments due to the different specifications of prior flux parameter error correlations; the error bars on the prior fluxes depict the average prior uncertainties. The reductions in the posterior estimated uncertainties in domain-integrated NEE relative to the prior uncertainties are 82.1% for the ALT_ECO experiment, 73.1% for ALT_DIST, and 77.6% for ALT_HYBRID. For the ALT_DIST and ALT_HYBRID experiments, the relative uncertainty reductions are a bit lower but generally in agreement with the relative error reductions. The ALT_ECO inversion results tend to underestimate the true errors, as was also seen when comparing Figures 7g and 7j. Both the prior and posterior uncertainties in all experiments underestimate the errors in regional NEE in the S and SE regions. In the other regions, as well as for the domain-wide NEE, the posterior uncertainties from the ALT_DIST and ALT_HYBRID inversions correspond well to the true errors, and the true values fall within two standard deviations of the posterior estimates.

Summary and Conclusions
This study introduced the TRACE Regional Atmosphere-Carbon Ensemble (TRACE) data assimilation system and demonstrated its CO 2 flux inversion capabilities in a series of controlled experiments. The main novelty of TRACE is its ensemble-based simultaneous state and parameter estimation (ESSPE) approach, which opens up possibilities to assimilate dense high-resolution observations from for example, spaceborne CO 2 observations. Here we focused on rigorously testing the system using synthetic observations from the conventional CO 2 tower network. Some error sources were deliberately excluded to simplify the interpretation of the results. As a consequence, the results are overly optimistic, and should be considered as a demonstration of the capabilities of the TRACE framework, and not as an assessment of the observational network.
The experiments were performed for July 2016 in a domain covering a large part of North America at a spatial resolution of 27 km. Most inversion experiments assimilated observations using a 1-week assimilation window. Under ideal conditions where the true flux error covariances are known, TRACE was able to reconstruct the true fluxes to a large degree; the posterior fluxes reduced the errors in monthly domain-integrated NEE by about 97% relative to the prior flux errors, and with spatial correlations of r = 0.87 and above between the posterior and true fluxes. In the second set of experiments, we considered a more realistic case where the true flux errors and their statistics are unknown. The true fluxes in these experiments were derived from different datasets than the prior fluxes, and three different prior flux parameter error correlations were tested. The inversions with distance-based and hybrid (ecoregion-and distance-based) flux error correlations resulted in the largest error reductions in monthly domain-integrated NEE (88.5% and 84.4% relative to the prior errors, respectively). In terms of regional fluxes on a spatial scale of 10 6 km 2 , the distance-based inversion was outperformed by the hybrid-based inversion (41.3% and 50.9% relative error reductions, respectively). The inversion with ecoregion-based flux parameters performed the worst for both domain-integrated (80.6% relative error reduction) and regional (34.5% relative error reduction) NEE, but nevertheless substantially reduced the biases in the prior fluxes.
TRACE provides not only estimates of CO 2 fluxes, but also uncertainties on the flux estimates derived from the ensemble spread. The inversion with ecoregion-based flux parameters had a slight tendency to underestimate the uncertainties due to the erroneous specification of prior flux parameter error correlations. For the distance-based and hybrid-based inversions, the estimated posterior uncertainties in domain-integrated and regional NEEs correspond well to the true errors, except in regions where the prior uncertainties were already underestimated.
In the ESSPE approach, the state vector consists of both CO 2 flux parameters and atmospheric CO 2 mole fractions. The CO 2 mole fractions are updated jointly with the flux parameters to reflect the adjustments to the parameter values. In this formulation, assuming persistent flux parameter errors, it is not necessary to include an assimilation window to optimize flux parameters that are distant from the observations, as the sensitivities of atmospheric CO 2 mole fractions to CO 2 flux parameters are carried by the ensemble perturbations and maintained between assimilation cycles. An assimilation window can be used to smooth the inversion solution over the time window, and to account for flux parameter errors that are decorrelated in time. A possible concern with including atmospheric CO 2 mole fractions in the state vector is that the mass of atmospheric CO 2 is not necessarily conserved, especially when covariance localization is applied. However, as shown by the experiments in this study, these errors were not detrimental to our regional inversion results, possibly because the errors in atmospheric CO 2 mole fractions were continuously transported out of the limited-area model domain, and thus did not accumulate over time (see e.g., H. W. Chen et al., 2019). For global dual-state inversions it may be beneficial to introduce algorithms that reduce mass conservation errors (Z. Liu et al., 2022).
This paper focused on giving an overview of the TRACE framework and ESSPE approach, and did not explore the full potential of the coupled atmosphere and carbon data assimilation system. Nevertheless, some new innovations were already shown, such as the dynamic selection of observations and hybrid ecoregion-and distance-based prior flux parameter error correlation scheme. Future studies will consider other important aspects of CO 2 inversions including assimilation of XCO 2 from satellites, more refined prior flux parameter error covariances, optimization of CO 2 lateral boundary conditions and fossil fuel emissions, inclusion of atmospheric transport errors, and assimilation of meteorological observations to improve the simulated atmospheric transport. Here we have shown that the ESSPE approach is suitable for regional CO 2 flux inversions, and that the TRACE system offers a promising platform for developing next-generation coupled atmosphere and carbon data assimilation techniques.

Data Availability Statement
All data used in this study are publicly available.