A Regional High-resolution Carbon Flux Inversion of North America for 2004

Resolving the discrepancies between NEE estimates based upon (1) ground studies and (2) atmospheric inversion results, demands increasingly sophisticated techniques. In this paper we present a high-resolution inversion based upon a regional meteorology model (RAMS) and an underlying biosphere (SiB3) model, both running on an identical 40 km grid over most of North America. Current operational systems like CarbonTracker as well as many previous global inversions including the Transcom suite of inversions have utilized inversion regions formed by collapsing biome-similar grid cells into larger aggregated regions. An extreme example of this might be where corrections to NEE imposed on forested regions on the east coast of the United States might be the same as that imposed on forests on the west coast of the United States while, in reality, there likely exist subtle differences in the two areas, both natural and anthropogenic. Our current inversion framework utilizes a combination of previously employed inversion techniques while allowing carbon flux corrections to be biome independent. Temporally and spatially high-resolution results utilizing biome-independent corrections provide insight into carbon dynamics in North America. In particular, we analyze hourly CO 2 mixing ratio data from a sparse network of eight towers in North America for 2004. A prior estimate of carbon fluxes due to Gross Primary Productivity (GPP) and Ecosystem Respiration (ER) is constructed from the SiB3 biosphere model on a 40 km grid. A combination of transport from the RAMS and the Parameterized Chemical Transport Model (PCTM) models is used to forge a connection between up-wind biosphere fluxes and downwind observed CO 2 mixing ratio data. A Kalman filter procedure is used to estimate weekly corrections to biosphere fluxes based upon observed CO 2. RMSE-weighted annual NEE estimates, over an ensemble of potential inversion parameter sets, show a mean estimate 0.57 Pg/yr sink in North America. We perform the inversion with two independently derived boundary inflow conditions and calculate jackknife-based statistics to test the robustness of the model results. We then compare final results to estimates obtained from the CarbonTracker inversion system and at the Southern Great Plains flux site. Results are promising, showing the ability to correct carbon fluxes from the biosphere models over annual and seasonal time scales, as well as over the different GPP and ER components. Additionally , the correlation of an estimated sink of carbon in the South Central United States with regional anomalously high precipitation in an area …


Abstract.
Resolving the discrepancies between NEE estimates based upon (1) ground studies and (2) atmospheric inversion results, demands increasingly sophisticated techniques. In this paper we present a high-resolution inversion based upon a regional meteorology model (RAMS) and an underlying biosphere (SiB3) model, both running on an identical 40 km grid over most of North America. Current operational systems like CarbonTracker as well as many previous global inversions including the Transcom suite of inversions have utilized inversion regions formed by collapsing biome-similar grid cells into larger aggregated regions. An extreme example of this might be where corrections to NEE imposed on forested regions on the east coast of the United States might be the same as that imposed on forests on the west coast of the United States while, in reality, there likely exist subtle differences in the two areas, both natural and anthropogenic. Our current inversion framework utilizes a combination of previously employed inversion techniques while allowing carbon flux corrections to be biome independent. Temporally and spatially high-resolution results utilizing biome-independent corrections provide insight into carbon dynamics in North America. In particular, we analyze hourly CO 2 mixing ratio data from a sparse network of eight towers in North America for 2004. A prior estimate of carbon fluxes due to Gross Primary Productivity (GPP) and Ecosystem Respiration (ER) is constructed from the SiB3 biosphere model on a 40 km grid. A combination of transport from Correspondence to: A. E. Schuh (aschuh@atmos.colostate.edu) the RAMS and the Parameterized Chemical Transport Model (PCTM) models is used to forge a connection between upwind biosphere fluxes and downwind observed CO 2 mixing ratio data. A Kalman filter procedure is used to estimate weekly corrections to biosphere fluxes based upon observed CO 2 . RMSE-weighted annual NEE estimates, over an ensemble of potential inversion parameter sets, show a mean estimate 0.57 Pg/yr sink in North America. We perform the inversion with two independently derived boundary inflow conditions and calculate jackknife-based statistics to test the robustness of the model results. We then compare final results to estimates obtained from the CarbonTracker inversion system and at the Southern Great Plains flux site. Results are promising, showing the ability to correct carbon fluxes from the biosphere models over annual and seasonal time scales, as well as over the different GPP and ER components. Additionally, the correlation of an estimated sink of carbon in the South Central United States with regional anomalously high precipitation in an area of managed agricultural and forest lands provides interesting hypotheses for future work.

Introduction
Carbon dioxide inversion studies have generally been focused on improved estimation of terrestrial carbon fluxes such as Ecosystem Respiration (ER), Gross Primary Production (GPP), and Net Ecosystem Exchange (NEE) as a means to better understand the carbon cycle of the earth. Researchers have progressively increased the resolution, in both time and space, and accuracy of the carbon flux estimates Published by Copernicus Publications on behalf of the European Geosciences Union.

1626
A. E. Schuh et al.: A regional high-resolution carbon flux inversion over the past decade. Early inversion studies were focused primarily with finding an explanation for the missing sink of carbon that can be easily identified from calculating a budget from annual fossil fuel emissions to the atmosphere, the effect of land use changes, and the oceanic carbon sink and comparing it to annual records of increasing atmospheric carbon dioxide concentrations. Given that the sink often represents a third of the annual fossil fuel emissions, it is of great interest to scientists and policy makers alike. Inversion results have been very effective at identifying large defining features of the terrestrial portion of the carbon sink (Fan et al., 1998;Gurney et al., 2002) although much debate remains even at extremely large scales (Stephens et al., 2007). However, the debate on a global scale has not deterred researchers from focusing these techniques on finer scale problems. In fact, criticism has been aimed at large scale global inversions because of the fact that their estimates can be biased on finer regional scales (Kaminski et al., 2001). The data available for regional inversion studies is increasing rapidly year after year, primarily within the developed industrial nations of the Northern Hemisphere. This provides researchers with some of the first opportunities to perform inversion studies in a very high-resolution setting. Gerbig et al. (2003) provided the first major regional inversion paper. They used a receptor-oriented inversion approach to investigate a series of flights from the CO 2 Budget and Rectification Airborne (COBRA) study conducted in 2000. Results showed that the effect of biosphere carbon fluxes could be seen at altitude in mixed layer CO 2 observed by aircraft. The paper pointed out several areas for future improvements in regional inverse modeling including improving biosphere-atmosphere exchange and convective transport modeling. Peylin et al. (2005) followed this with a regional inversion based on western Europe in which he estimated daily fluxes for a month using relatively continuous measurements of CO 2 from towers in the inversion domain. The most similar effort made for North America comes from the ongoing CarbonTracker project . Peters et al. (2007) used a nested transport structure (TM5) with a relatively high-resolution 1-degree inner grid over North America. A priori carbon fluxes were estimated by modifying 1-degree by 1-degree monthly output from the Carnegie Ames Stanford Approach (CASA) model to provide diurnal variability by incorporating a Q 10 temperature relationship for respiration and a linear scaling of GPP with solar radiance. NEE estimates were optimized by estimating linear correction factors for NEE for each of up to 19 ecoregionbased (Olsen et al., 1992) sub-areas of North America based upon a 5-week smoothing window. The coarseness of the inversion over North America is required in order to regularize the inversion problem in light of limited observations.
Our inversion framework has drawn upon certain techniques from previous inversions while including some new features. The aim of the inversion is to provide fine scale inversion results over North America for 2004. A novel feature of this inversion is the distinct estimation of GPP and ER instead of just NEE, which to our knowledge has not previously been performed, at least in the regional framework. We have drawn upon the spatial correlation constraints used by Rödenbeck et al. (2003) and Michalak et al. (2004), largely in order to regularize the inversion problem. We are largely using the spatial correlation in the prior error covariance structure to regularize the problem. Attempts were made to estimate the spatial correlation via the measurements but the data was not constraining enough. Future alternatives that might be possible would be to estimate parameters via a prior flux model as in Mueller et al. (2008) or . However, it has been shown in practice that certain isotropic spatial error correlations can work well as a regularization agent. Carouge et al. (2008a, b) investigated the potential of a 10 tower network of CO 2 observing towers over Europe using a 50 km resolution grid over Europe. They found that 10 days temporal and 1000 km spatial averaging was required in order to obtain good agreement between estimate and "true" fluxes. Surprisingly, they found that these isotropic assumptions on the spatial errors performed better than an estimate of the spatial errors based on the "physical errors", those that could be calculated by knowing the "true" fluxes.
Large matrix inversions limited the inversion grid resolution to approximately 10 000 km 2 (60⇥36 grid composed of 100 km by 100 km grid cells). For sensitivity studies involving numerous inversion runs, a 40 000 km 2 grid (30⇥18 grid composed of 200 km by 200 km grid cells) is used. Many previous global inversions have been performed upon grid areas of around 5 to 10 times that size. In order to provide some contrast, CarbonTracker optimizes 17 bias correction factors for NEE over North America (with 4 of those representing less than 0.5% of the land area each) while this inversion typically optimizes 30⇥18=540 each (30⇥18 grid mentioned above) for ER and GPP. It is important to note that the employment of a spatial correlation constraint and decorrelation length scale, either due to a formal statistical model, or as a method of regularization, does reduce the effective degrees of freedom so that we certainly do not expect the optimization of 540 "independent" parameters. Nevertheless, it is important to note that Schuh et al. (2009) showed how biases can occur when using a relatively coarse fixed set of regions within an atmospheric inversion as opposed to a finer set of regions, even when assuming spatial-scale patterns of carbon flux on the order of 500 km and greater. The flexibility we have achieved by avoiding fixed inversion regions does not come without a cost since we cannot simultaneously optimize fluxes outside of North America. Therefore we used offline-derived boundary conditions and provided these as fixed contributions to the tower CO 2 budget. Schuh et al. (2009) showed that considerable success could be achieved in estimating large spatial scale ER and GPP signals in the midst of small spatial scale variability in fluxes. We leveraged this result and put the problem in a Kalman filter framework in order to allow higher resolution spatial estimation. This filter is of a somewhat simple variety and allowed us to work with all portions of the inversion, such as complete prior and posterior covariance matrices, explicitly. We then tested sensitivity to a number of pieces of the inversion considered uncertain, including parameters in the actual inversion as well as fixed contributions to the modeled CO 2 such as fossil fuel and boundary inflow. As far as we know this is also the first paper providing a comparison of inversion results derived by using two independent boundary inflow estimates. Additionally, the effect of including recently available high-resolution fossil fuel inventory data is quantified.

Prior flux model and transport
The Simple Biosphere model (SiB) is based on a land-surface parameterization scheme originally used to compute biophysical exchanges in climate models (Sellers et al., 1986), but later adapted to include ecosystem metabolism (Sellers et al., 1996a;Denning et al., 1996a). SiB has been coupled to the Brazilian version of the Regional Atmospheric Modeling System (RAMS; Pielke et al., 1992;Frietas et al., 2006) and used to study PBL-scale interactions among carbon fluxes, turbulence, and CO 2 mixing ratio  and regional-scale controls on CO 2 variations (Nicholls et al., 2004;Wang et al., 2006). This latest version of SiB is termed SiB3.
In SiB3, Net Ecosystem Exchange (NEE) is composed of two component fluxes, Gross Primary Productivity (GPP) and Ecosystem Respiration (ER), which includes autotrophic (canopy respiration and root respiration) and heterotrophic respiration terms (due to decomposition of dead organic matter), NEE(x,y,t) = ER(x,y,t) GPP (x,y,t) (1) where x and y represent grid coordinates and t represents time. High-frequency time variations of photosynthesis and respiration are assumed to be well understood and easily modeled processes, i.e. due to diurnally varying quantities such as radiation, temperature, or longer term variations in modeled quantities such as soil moisture etc. Photosynthesis and assimilation are derived using a coupling of equations based upon the work of Farquhar, Collatz, and Ball (Farquhar et al., 1980;Collatz et al., 1992;Ball et al., 1987) while soil respiration is based upon a rather simple function of temperature and soil moisture and constrained in such a way that annual NEE is equal to zero (Raich et al., 1991;Denning et al., 1996). Several papers have provided comparisons of models to observations, largely by using eddy flux towers to measure true fluxes of water, carbon, and energy (Baker et al., 2003, Water ( 2008; Hanan et al., 2005). Longer-term, more persistent biases are estimated by solving for unknown multiplicative biases in each component flux after smoothing in space and time. While these biases could result from incorrectly modeled short term processes, such as errors in the daily development of the planetary boundary layer, or short-term processes not in the model such as seasonal fertilization and irrigation, the main purpose is to capture longer-term processes not explicitly modeled such as land use change (Robertson et al., 2000;Peterson et al., 1998), disturbances, anthropogenic fertilization effects (Oren et al., 2001), managed forestry (Tillman et al., 2000), and large scale carbon removal . This modeling is accomplished by convolving the influence functions generated from a lagrangian particle dispersion model, LPDM (Uliasz and Pielke, 1991;Uliasz, 1993Uliasz, , 1994Uliasz et al., 1996;Zupanski, 2007), with gridded Gross Primary Productivity (GPP) and total Ecosystem Respiration (ER) at each time step in SiB3-RAMS. The LPDM transport scheme reverses advection derived from RAMS at very fine time scales and parameterizes vertical turbulent diffusion according to a Gaussian process. A large advantage of this model is the ability to simulate transport of atmospheric constituents at sub grid scales, reducing representation error that might be caused by associating an observing tower with a 40 km grid cell in the model. By tracking particles upwind, backward in time, from the towers, one may make inferences about the contribution of upstream GPP and ER sources.
In particular, we have estimated regional fluxes from atmospheric mixing ratios by assuming that the model of the component fluxes is biased, and that the biases are smoother in time and space than the fluxes themselves: NEE(x,y,t) = (1 + RESP (x,y)) ER(x,y,t) (1 + GPP (x,y)) GPP(x,y,t) The model domain, shown in Fig. 1, consists of most of the United States as well as a large portion of Canada and the northern portions of Mexico. Both SiB3 and RAMS were run on a single 150⇥90 grid of 40 km cells, with SiB3 utilizing 3 patches per cell to capture subgrid-scale variability in land cover. RAMS meteorology was nudged with 40 km forecast meteorology from the National Center for Environmental Protection's Eta model throughout the domain using a 4 Dimensional Data Assimilation (4DDA) scheme to produce more reliable wind fields. Soil classes were calculated from 5 min "%clay/% sand/% silt" soil data from the International Geosphere-Biosphere Programme (IGBP). Biomes were extracted from the UMD classification scheme of the MODIS 12 Landcover 1 km product and mapped to the most similar SiB biome class for all cells and for each of the three patches used. An exception are the C 4 vegetation classes, grasses and crops, which were projected onto the MODIS biomes from (Wang et al., 2006). The crop characterization is admittedly simple and more work is currently being done to incorporate more accurate crop maps and more realistic crop modeling into SiB (Lokupitiya et al., 2009). SiB has traditionally calculated fPAR, which defines the fraction of photosynthetically available radiation that is absorbed by the plant canopy, and Leaf Area Index (LAI) using satellite derived NDVI fields. The code was changed to use fPAR and LAI fields derived by the Moderate Resolution Imaging Spectroradiometer (MODIS) (Mu et al., 2007) and averaged over appropriate biome-areas based upon the three patch scheme. SiB3 was run with these 8-day fPAR and LAI products that were provided by the Numerical Terradynamics Simulation Group at the University of Montana who generated it for use in constructing the official Moderate Resolution Imaging Spectroradiometer GPP product.
Modeled carbon dioxide at the tower is calculated as the sum of 3 component fluxes convoluted by time and tower dependent transport. The boundary inflow component was calculated by convolving the influence functions from the LPDM model over boundary CO 2 fields derived using a global biospheretransport model. At any point in time, the boundary inflow is the average of all upstream particles located in a 3 dimensional 40 km thick rectangular "ring" around the domain. CO 2 resulting from the transport of fossil fuels to the towers is calculated by convolving the influence functions from the LPDM model with surface fossil fuel flux estimates. In particular, the boundary CO 2 fields were calculated by combining transport from the Parameterized Chemistry Transport Model (PCTM) (Kawa et al., 2004;Parazoo, 2007) and precalculated archived hourly SiB3 fluxes (Baker et al., 2007) on a 1.25-degree by 1-degree global grid. The model was spun up for 2000-2004 and the CO 2 was centered around the Northern Hemispheric mean CO 2 for 2004. In addition to this, results from the CarbonTracker project, which provide globally optimized CO 2 concentration fields, are used for comparison purposes.
Fossil fuel fields were constructed using recently available high resolution Vulcan fossil fuel inventory fields (Gurney et al., 2008), at a 10 km horizontal spatial scale and hourly temporal scale. Previously available fossil fuel flux fields were derived by distributing country-level fossil fuel sources spatially as a function of population at a 1-degree resolution (Andres et al., 1995). The Vulcan fields provide many improvements including the incorporation of mobile emission sources and power plants, often located in areas distant from high density population centers, increased temporal resolution allowing the modeling of diurnal variability, and increased spatial resolution allowing better delineation of high density population centers. The sensitivity to the new fossil fuel fields is tested by running inversions using both the Vulcan fields as well as the Andres et al. (1995) fields.
The effect of this on boundary inflow estimates is that the PCTM-SiB3 calculated boundary CO 2 fields lacks the effect of sources or sinks in 2004. Given the consensus opinion of an annual mean sink for carbon resulting from the biosphere, this means that the CO 2 fields used will be biased somewhat by the effect of not including this expected global sink. We investigate the effect of this by including a comparison of the inversion using CarbonTracker optimized CO 2 concentration fields for boundary inflow, which includes the effects estimated sources/sinks outside of the regional modeling domain. As of this time, carbon dioxide resulting from forest fires is not included in the global PCTM-SiB3 inflow or domain SiB3 runs, but is included in the CarbonTracker inflow providing one more contrast between the two fields.

Observational data
Half-hourly averaged calibrated CO 2 observations were provided for eight measuring sites (WLEF, Harvard Forest;Urbanski et al., 2007; ARM, BERMS, Fraserdale, Western Peatland, WKWT, and Argyle -ME) for 2004 (Parazoo, 2007). Gerbig et al. (2003) found mean standard deviations on the order of 0.6 to 1 ppm when viewing morning and afternoon vertical profiles of CO 2 in the mixed layer. As a consequence, robust afternoon snapshot observations, at 12, 2, 4, and 6 p.m. LT, are used in order to avoid inversion model sensitivity to poor atmospheric transport modeling of extremely stable and stratified nocturnal atmospheric conditions near the ground. One exception is the WKWT tower in Moody, TX. For most days, data at this tower consistently showed high CO 2 concentrations in the 12 p.m. LT records that were more consistent with typical morning CO 2 than with wellmixed afternoon CO 2 . For this tower, mixed boundary layer conditions appeared to be better represented by snapshot observations shifted by 2 h: 2, 4, 6, and 8 p.m. LT. The first 10 days of the year are not comparable due to a lack of transport preceding 2004. In all there were 2433 missing observations, resulting in 4 (observations/day) ⇥ 8 (towers) ⇥ 355 (days) 2433 (missing) = 8927 observations.
In a previous pseudo-data inversion using a very similar model (Zupanski et al., 2007), the errors on the observations were assumed to be 1 ppm for afternoon observations. Nevertheless, relative to the inversion techniques presented in the next section, the errors on these observations should include errors due to calibration error, mapping error, transport error, and representation error. For this inversion, transport error and representation error are likely the largest components which are notoriously tricky to quantify. Investigations into the sensitivity of inversion test results combined with initial maximum likelihood estimation results suggest errors in the range of 5-6 ppm are appropriate for this particular inversion. For the remaining inversions, the errors are assumed to be identical and independently distributed (i.i.d.) mean zero errors with standard deviation set to 5.5 ppm. This is a simple assumption and we certainly do not expect to the error to be completely homogeneous across towers although at what scale the observation error should be estimated is still somewhat uncertain. We also note that there certainly is expected to be autocorrelation in the errors within a daily time frame so that the "effective" number of observations is likely much less than 4 each day. The end result is that the observational error term over multiple observations is probably estimated as being somewhat lower than reality. For example, a mean of 4 afternoon observations has an estimated 2.75 ppm error, based on Gaussian 5.5 ppm independent errors for each observation. In reality, the error of the mean observation is probably larger due to likely temporal correlation in the observation errors.

Climatic conditions for 2004
The 2004 year was the 6th wettest in the contiguous United States over the preceeding 110 years . It was also warmer than on average. Nevertheless, there was a great amount of variability in precipitation and temperature as a function of location and season. Drought continued in the west through the summer of 2004, essentially prolonging a multi-year period of drought conditions. The spring was also very dry for the southeast, extending a period of dry conditions from late in 2003. However, summer brought increased precipitation to the east and southeast, culminating in enormous amounts of rain in late summer and early fall due to an extremely active hurricane season. The south (Texas, Louisiana, Mississippi, Arkansas, Oklahoma, and Kansas) had the wettest summer on record and was much cooler than average. These conditions were important as they provided initial conditions for SiB3 that involved soil moisture induced plant stress over large areas of the United States.

Inversion technique
Standard multivariate Gaussian assumptions are made and data are assimilated using a modified Kalman Filter algorithm (Kalman, 1960). In particular, for an initial length n CO 2 measurement vector y representing the first set of measurements, length m unknown CO 2 flux bias vector (dimensionless), n⇥n observation error covariance matrix 6 (ppm 2 ), n⇥m Jacobian flux-transport matrix G (ppm), length m prior flux bias estimate 0 (dimensionless), and m⇥m model-prior mismatch covariance matrix 6 0 (dimensionless), the Bayesian statistical assumptions are a Gaussian distribution on the "measurement" errors as well as a Gaussian distribution on the a priori distribution of , i.e.: The posterior distribution of the flux bias vector can be solved for analytically and is: With a little bit of algebra, one can rewrite the mean of the posterior distribution of the mean, giving the Kalman-filter updating equation for the mean.
The posterior mean and variance of x are then fed into the next filter step with a new set of measurements. This particular inversion estimates biases over 7-day periods using available data from that 7-day period of time. Therefore, bias estimates for both ecosystem respiration and GPP as well as corresponding variance estimates are available for all of 2004 with the bias estimates changing with a weekly resolution.
Two difficulties often arise when using filter-style correction schemes. The filter estimates can drift away from realistic values if the data are not plentiful or precise enough to constrain it. Secondly, the nature of the Kalman filter at each step is to create posterior variance estimates that are in general smaller than the prior estimates. This can essentially cause the filter to get "stuck", when an explicit dynamical model of the biases is not available, and thus produce unrealistically small posterior variance estimates around the biases. There is generally no easy solution to this problem. Artificially inflating the posterior variance at each filter step is one method in which one can try to circumvent (Zupanski et al., 2007). This accommodates the fact the biases are likely to 1630 A. E. Schuh et al.: A regional high-resolution carbon flux inversion change in reality and it allows the filter to consider a wider range of possibilities for the bias factors. However, it does not necessarily constrain the biases to any particular "reasonable" region of values allowing the bias estimates to drift into unrealistic parameter space. Therefore, we have chosen to weight the filter at each step with a "grand" prior. This effectively handles both of the preceding problems. With respect to our inversion, there will be three pieces of information at each step, the grand prior which is derived from the forward SiB3-RAMS model with an error assumption, the local prior which is derived from the previous filter step's posterior flux bias distribution, and the data which forms the statistical likelihood function. In some sense, this new piece of the covariance structure provides a bound upon how much the inversion can "learn" about the bias structure.
In order to quantify, we denote the grand prior as a multivariate Gaussian distribution around grand with covariance matrix 2 grand 6 grand , and additional weight factor w, and we rewrite the expression given in Eq. (4) as: Thus is distributed as a multivariate Gaussian with parameters: Equation (7) specifically separates out the variance scalars, 2 grand , 2 0 , and 2 obs from the covariance matrices, leaving the covariance matrices essentially scaled to 1. The w weight is a redundant factor and is simply included to facilitate easier interpretation of tightening/loosening of the grand prior covariance (around the SiB3 derived a priori carbon fluxes). Unless otherwise specified, this weight, w, on the grand covariance matrix is set to 2. This means that the initial variance around the grand prior is increased, thus providing a weaker constraint. For the initial filter step, only the grand prior is used. After that point, there exist both a grand prior and a prior (from the posterior of the previous filter step). The inversion is further constrained by the assumption of spatially correlated errors in the grand prior, i.e. the covariance matrix 6 grand will take on the following form.
The respiration and GPP covariance matrices are each formed from the exponential covariance function, where t i,j is the distance between points i and j . Cov The h 0 parameter is the decorrelation length scale parameter, giving the distance at which the covariance between two points is equal to 2 0 (1 ↵ 0 )e 1 . The 2 parameter is the scalar variance parameter and determines the variance of the marginal distribution of the particular flux component. The parameter ↵ 0 controls what percentage of the covariance can be attributed to spatial covariance, as opposed to spatially independent errors, often termed "nugget" variance. While the "nugget" parameter is an important parameter if one is fitting a rigorous statistical spatial model to the errors, for regularization purposes it is often set to zero which is what we will do for the remainder of the paper.
It is important to note that the use of a high resolution grid for the inversion certainly does not imply that meaningful inferences can be made at the finest scale. Assumed decorrelation length scales of 500 km and greater certainly imply a strong constraint upon the solution effectively giving somewhat smoothed solutions spatially. For example, while 540 parameters are used spatially, there are only 8 towers providing information at any point in time (assuming no missing data) and those are only afternoon observations. In effect, if the wind is coming from N.W. portion of the domain to all the towers, then the inversion can only learn about the N. W. portion of the domain and only within the confines of the differences in upstream sampling footprints and the differences in the observed CO 2 at the towers. For example, based upon a 1000 km decorrelation length scale smoothing scheme with a two week assimilation cycle, the effective degrees of freedom of the data used in the paper might only be between 2 and 7, with the E dimension being estimate between 15 and 30 (Park and Xu, 2009). This can be visualized by noticing the scale of the corrections in many of the figures in the paper.
It was shown in Schuh et al. (2009), that under isotropic error type conditions (for assumed and "true" errors) that this inversion model is robust to small spatial scale random deviations in flux bias and that post-aggregated (in space) estimates can be very good even when using a fairly sparse network of towers observing CO 2 . The usefulness of the high resolution grid allows one to separate the impact of a tower residual across space more effectively. For example, if a tower only saw a corner of a somewhat large inversion region than the correction imposed on that corner is essentially projected onto the entire inversion region. This essentially arises from a lack of uniform sampling over every inversion "grid cell". The effect of this potential error is of course dependent upon the "actual" structure of the errors which is not often known. Nevertheless, regardless of whether one chooses to assume isotropic error conditions through a spatial smoothing such as ours, or larger inversion regions similar to Peters et al. (2007), given the unconstrained nature of the inversion problem, it is always important to assess the impact of varying certain unknown parameters in the inversion, such as the choice of inversion regions: grid vs biome, localization schemes, the spatial decorrelation length scales, the weight given to the "grand" prior, and the fixed CO 2 contributions from both the boundary inflow and fossil fuel sources.

Sensitivity
The inversion essentially guarantees some improvement in prediction of observed CO 2 (Eq. 5). However, when using a regression style approach in a heavily unconstrained environment, this improvement can often be overstated because of the great freedom the inversion has to fit the data. Therefore, it is often desirable to go beyond simply comparing observed carbon dioxide at the towers to model-based predicted carbon dioxide. Comparing model observations to independent observations not used in the inversion, comparing models which predict similar quantities, as well as testing the sensitivity of the model to variations in unknown parameters are all methods of generating more confidence in estimates.
We used a variety of different procedures to test the sensitivity of the inversion. Therefore, we first test the sensitivity of the inversion to varying the inflow of CO 2 at the boundaries. To do this, we derive boundary inflow to the 8 towers using the LPDM model and optimized carbon dioxide concentration fields from the CarbonTracker project . Inversion results are then compared with the results derived from the LPDM model and the PCTM inflow. Secondly, we vary several different variance parameters and derive annual domain-summed NEE and tower observation based RMSE based upon the varied parameters. Thirdly, we use a re-sampling procedure in which we create 45 different observation data subsets by holding out a randomly selected 50% of the observation data for each. Each set of data is run through the weekly inversion scheme and the sensitivity of the predicted CO 2 at the towers and the estimated flux biases is explored. This provides estimates of the variability of the flux correction factors and can be used to assess the sensitivity of the source/sink to the constraint provided by the data. Using the held out data as independent evaluation data and the complementing data as training data for the inversion, one may also derive a more accurate estimate of Root Mean-Squared Error (RMSE) of the inversionoptimized fluxes. We test the impact of the high resolution Vulcan fossil fuel inventory on the inversion results by comparing inversion results relying upon Vulcan to those results utilizing the Andres et al. (1995) fossil fuel inventory.
SiB3 has been evaluated at many sites and over many time periods, nevertheless, the particular model run used for the a priori flux estimates was not optimized to fit the flux data at any site in particular. Even though there is a mismatch in representation, with the flux towers representing footprints of less than a square kilometer and the inversion results representing flux estimates on the scale of thousands of square kilometers, we believe that these comparisons are of value, especially in locations that are more spatially homogeneous than others, such as grasslands and large forest reaches. This is then the fourth comparison we make.

Results
As was indicated in the previous section, there are a number of variables that the inversion will likely be sensitive to and therefore the results are expected to be quite variable. For results, we choose to present one particular case with a fixed set of inversion inputs as an initial case study and then use it to compare the effect of varying the boundary inflow and the source of the domain fossil fuel fluxes. With reference to the preceding section and Eq. (7) in particular, the following values are used for these inversions: grand =0.25, 0 =0.25, obs =5.5 ppm, w=2, h 0 =1000 km. In particular, a value of grand =0.25 would mean that we expect that approximately 68% of the GPP and ER biases are within ±25% of the original SiB3 estimated fluxes, with 95% within ±50%. This variation when combined with positive spatial correlations was shown to provide a reasonable a priori range of annual domain-summed NEE. These deviations must generally be kept to less than 30-40% to ensure that posterior ER and GPP fluxes are not reduced by more than 100%, which makes no conceptual sense. We then test the sensitivity of the results over a number of varying inversion inputs using the PCTM boundary conditions and the Vulcan fossil fuel flux field.

General structure of results
CO 2 can be predicted by invoking the relationship shown in Eq. (3). The predicted mean observed CO 2 is derived as Gx wherex represents one (for the prior fluxes) plus the inversion-optimized flux biases. Using the PCTM boundary conditions and the Vulcan fossil fuel inventory, a comparison of the inversion-corrected posterior predictions at the towers to the observations is shown in Fig. 2. For domain-summed temporal plots, NEE is calculated via Eq. (2) while ER and GPP are calculated via the two respective summands on right hand side of that equation. These domain-summed temporal results are shown in Fig. 3.
The observed carbon dioxide concentrations contain information that infers a dampening of the a priori annual GPP cycle, and hence the a priori annual ER cycle (due to the strong correlation of the annual sums of each). Since both GPP and ER are significantly dampened, it is not surprising that the NEE signal is dampened as well. Furthermore, the data suggest a weak temporal shift in the prior NEE signal. This manifests itself as a stronger, but more gradual onset     of spring, followed by a weaker overall carbon sink over the middle and late summer periods. It was brought to our attention by Steve Wofsy (Harvard U.) that this discrepancy in a priori and a posteriori GPP and respiration fluxes could arise from a bias in the meteorological data that was driving SiB3 to produce the a priori fluxes, in particular biases in shortwave radiation. Upon investigation and comparison of several different reanalysis products to recently available Ameriflux radiation observations, it does appear that nearly all of the reanalysis products investigated had somewhat uniformly, and far from insignificant, positive biases in shortwave radiation. This certainly could play a role in the a priori GPP and respiration fluxes being significantly larger than they should be. This is currently being investigated under the guise of the North American Carbon Program (NACP) by Daniel Riccuito (Oak Ridge National Laboratory), and a manuscript is currently in preparation describing further the differences found between observations and reanalysis products and the possible effect on biospheric carbon flux estimates. Nevertheless, this is an a priori estimate of fluxes, and although we certainly would like it to be close to the posterior estimate, it simply represents the best knowledge we currently have about the processes.
We use the resampling procedure, that was first mentioned in Sect. 2.5, to account for variability that might be associated with over fitting the model and which provides additional variability to the standard covariance estimates of the biases given in Eq. (6). Forty-five different inversions are run, each based upon a different subsample of the observations. Assuming temporal independence of the errors in the filter, one may simulate properties of the annual NEE probability density functions (pdf) for each of these 45 inversions by using the posterior covariance provided at each step of the Kalman Filter for each inversion. A 95% Confidence Interval (CI) for the entire domain can be calculated at each step of the filter for each of the 45 inversions. The CI shown in Fig. 3 then characterizes variability in the NEE by selecting the 95% CI of each set of 95% CIs for each weekly time step.
The ensemble mean of the domain summed annual NEE flux is approximately 0.68 Pg/yr while the standard deviation of this estimate is about 0.11 Pg/yr. It is important to note that this standard deviation estimate does appear to be too small, giving tighter bounds on the flux than found in other inversion papers Peters et al., 2007). An additional source of variability in the estimate is discussed later (Sect. 3.4) and likely provides another 0.1-0.15 Pg/yr to this standard deviation estimate. The spatial representation of these sources and sinks can be seen in the first panel of Fig. 7. Depictions of this variability in a spatial framework are shown in Fig. 4. This variability is partitioned into two pieces, variability associated with the spread of mean estimates over the 45 inversions (measure of over fitting) and variability associated with summing up the posterior variances at each filter step (regular KF variance)  Finally, the square root of this summed variance (standard deviation) is displayed and is a measure of the uncertainty of the mean estimate due to model over-fitting. For the right panel, the summed annual variance in NEE is calculated for each inversion, from the weekly filter estimates, and the the square root of this (standard deviation) is shown for each cell. These plots aim to provide a measure of the uncertainty of each cell's NEE estimate, incorporating the correlation between ER and GPP in each cell, but not incorporating the spatial correlation in the covariance matrices. evaluated over all 45 inversions. Besides the spatial display of posterior variance information for NEE, which roughly tracks the convolution of the sampling footprint of the network and the prior ER/GPP signals, the results show that over fitting the model may provide a significant source of variability comparable to that which is normally constructed from each filter step's posterior covariance matrix.
The residuals from the model fit are generally symmetric and do not appear to deviate substantially from normality (Fig. 5). There is a slight but pronounced positive skew to the residuals indicating that when the residuals deviate most strongly from zero, the observed CO 2 is greater than the modeled CO 2 . Biases remain in the inversion process, likely a result of residual pdfs' deviations from symmetry. Sites in the north and northeastern portion of the domain appear to be most sensitive to this, in particular AMT and Harvard Forest (Table 1). These sites seem to be affected more by a strong a priori seasonal cycle than the other sites. Additionally, we note that the these towers are in relatively close proximity to the most populated areas of North America and it is possible that occasional spikes in anthropogenic emissions from the northeast coast of the United States could impact tower concentrations. Weekly chi-square statistics were calculated to diagnose the model's performance. Values near to one indicate that the assumed errors are being estimated reasonably, a priori. The weekly chi-square innovation statistics are generally near 0.5 from January through May and then around 1 for the summer and remainder of year. The innovation statistics show more temporal variability in the summer time. The low value in the winter time in indicative of some heterogeneity in the model-data residuals, seasonally, and that the assumed errors in the winter might be too large. In essence, this might weight the model too heavily towards the prior which might imply that magnitude of the winter time NEE adjustment, which is generally a C sink under our a priori flux scenario, might be too weak.

Sensitivity and robustness of results to inflow
Inflow of CO 2 from the boundaries has typically been a large concern of regional models (Gerbig et al., 2003;Peylin et al., 2005). In extremely limited domain problems, the variance of the CO 2 coming in from the boundary can easily dwarf the changes inside the domain due to local biotic uptake and release. Therefore it is of interest to gauge the sensitivity of the inversion to varying boundary inflows. The boundary conditions included in this model were constructed from a global simulation using SiB3 and PCTM (Parazoo et al., 2007). The CarbonTracker project has provided CO 2 mixing ratio data based upon globally optimized fluxes . SiB3 has no annual source/sinks whereas Car-bonTracker includes an annual source/sink estimated from observations of CO 2 . A plot of the difference between the two inflows is shown in Fig. 6. The inflow annual mean and temporal pattern is very similar for PCTM and Carbon-Tracker with the main difference being a seasonally stronger cycle in the PCTM-SiB3 results, likely a result of the underlying biosphere model, SiB3, providing a stronger seasonal GPP/NEE signal than the corresponding CASA model used in CarbonTracker. In addition to running comparison inversions between these two CO 2 inflow estimates, we also run the inversion with a fixed inflow estimate of 378 ppm representing the annually averaged PCTM inflow over the period of the simulation in order to show the necessity of reasonable boundary inflow values in calculating source/sink estimates. Figure 7 shows a comparison plot of maps of the annual mean NEE estimate based upon CarbonTracker (w/CASA), PCTM (w/SiB3), and the fixed inflow condition. The results are similar for the CarbonTracker and PCTM inflows. Both results have similar spatial and temporal characteristics but differ mainly in magnitude. The PCTM-based inversion results in a sink of 0.1-0.2 Pg/yr more than that of the CarbonTracker-based result. The PCTM-based boundary conditions do not account for the expected global carbon sink outside of the inversion domain, which forces the inversion to increase the North American sink to compensate. This results in the PCTM-inflow based inversion having a larger annual sink estimate than the CT-inflow based inversion. The sink estimated with the PCTM inflow was 0.65 Pg/yr while the sink estimated with the CarbonTracker inflow was estimated at 0.48 Pg/yr. It does seem somewhat surprising that the results from the two inflows are still close, within approximately 30% of one another. This indicates that local observations may be affected significantly more by local fluxes than by larger scale fluxes in distant locations outside of the model boundary.

Sensitivity of results to fossil fuel inventory
Until the release of the Vulcan fossil fuel inventory in 2008, most researchers were reliant upon the Andres et al. (1995) fossil fuel inventory, which was released at annual time 25 477

479
Inflow of CO 2 from the boundaries has typically been a large concern of regional models 480 (Gerbig et al., 2003;Peylin et al., 2005). In extremely limited domain problems, the 481 variance of the CO 2 coming in from the boundary can easily dwarf the changes inside the 482 domain due to local biotic uptake and release. Therefore it is of interest to gauge the 483 sensitivity of the inversion to varying boundary inflows. The boundary conditions 484 included in this model were constructed from a global simulation using SiB3 and PCTM 485 (Parazoo et al., 2007). The CarbonTracker project has provided CO 2 mixing ratio data 486 In particular, this figure shows the "difference" between estimates of CO 2 arriving at tower due to two distinct boundary inflows (1420 sequential '12/2/4/6PM' observation sequences for each of 8 towers.) ppm Fig. 6. Figure shows the effect of boundary inflow CO 2 upon tower CO 2 concentrations. In particular, this figure shows the "difference" between estimates of CO 2 arriving at tower due to two distinct boundary inflows (1420 sequential "12/2/4/6 p.m." observation sequences for each of 8 towers). scales and at a 1-degree resolution over the globe. For many large-scale inversion applications, this inventory is adequate. However, for higher resolution studies within the United States, the Vulcan fossil fuel inventory provides a dramatic improvement in both space and time accounting of fossil fuel fluxes. The main difference between these inventories is the redistribution of some fossil fuel sources from population centers to more distant locations representing mobile sources and power plants. The Vulcan fossil fuel flux estimates are at a much higher resolution in both time and space. Previous inversions had to grapple with the fact that some observing stations are located within enormous fossil fuel flux regions. For example, a semi-rural location like Harvard Forest would very likely be located in the same grid cell as the large metropolitan city of Boston. Given no sub-annual temporal resolution to the fossil fuel fluxes, an observing tower located at Harvard Forest was often seeing a 24 h continuous stream of fossil fuel fluxes arising from a city over 100 km away. However, the 10 km horizontal resolution of the Vulcan inventory allows these to be separated and additionally provides a diurnal and seasonal estimate of these fluxes, which is important for inversions based upon hourly observations.
In order to gauge the impact of incorporating the Vulcan data, we first contrasted the contributions to each of the 8 towers from each of the inventories. For many of the stations, the afternoon differences between the two were very small. Differences at the ARM site in Oklahoma, the WLEF site in Wisconsin, the Canadian sites, and the Argyle, Maine site were on the order of a few ppm. Differences at the Moody, Texas tower were in the range of 5 to 5 ppm. While the differences across most towers were relatively small, the differences at Harvard Forest were between 25 and 30 ppm! The difference in the annual NEE estimate is shown in Fig. 8. The effect on the inversion is far from trivial with differences of up to 150 g/m 2 per year recorded along the northeast coast of the United States, similar in magnitude to the maximum annual sinks estimated by the inversion. These differences are a result of coarse fossil fuel flux fields providing artificially high sources of CO 2 to the Harvard Forest tower which must be neutralized via a large local sink.

Sensitivity and robustness of results to prior variance structure
A test of the sensitivity and effect of the prior upon results is important because of the use of an informative Bayesian prior, that is, a prior flux estimate in which the inversion will likely be sensitive. With reference to Eqs. (5) and (7), the w, 2 0 , and h 0 parameters are varied and results are shown in Fig. 9. These figures show that results are sensitive to nearly all of these parameters, providing different degrees of RMSE and sink strength depending upon the particular combination. In particular, sink estimates range between 0 and 1 Pg/yr. The ensemble of estimates, over the various possible

566
A test of the sensitivity and effect of the prior upon results is important because of the use 567 of an informative Bayesian prior, that is, a prior flux estimate in which the inversion will 568 likely be sensitive. With reference to Eq. 5 and Eq. 7, the w, σ 0 2 , and h 0 parameters are 569 varied and results are shown in Fig. 9. These figures show that results are sensitive to 570 nearly all of these parameters, providing different degrees of RMSE and sink strength 571 depending upon the particular combination. In particular, sink estimates range between 0 572 and 1 Pg/yr. The ensemble of estimates, over the various possible a priori variance 573 parameters, has a standard deviation of approximately 0.2 PgC/yr. This likely contributes 574 another 0.1 PgC/yr to 0.15 PgC/yr (depending upon the existence of correlation between 575  Andres et al. [1995] fossil fuel inventory. Positive values (purple) indicate carbon sinks were stronger using Andres inventory. Spatially-summed annual difference between Vulcan-based NEE estimate for 2004 and Andres[1995] based NEE estimate for 2004 is less than 0.01 PgC. a priori variance parameters, has a standard deviation of approximately 0.2 PgC/yr. This likely contributes another 0.1 to 0.15 PgC/yr (depending upon the existence of correlation between the variance shown here and earlier variance estimates due to jackknife resampling and the Kalman filter posterior variances) to the initial standard deviation estimate of 0.11 Pg/yr given earlier. This would give an adjusted standard deviation estimate of approximately 0.2-0.25 PgC/yr to the posterior annual NEE estimate shown in Fig. 3.
An RMSE-weighted average of the sink estimates show a sink of 0.57 PgC/yr, 20% higher than our single case scenario that we have followed throughout these results. Values very near the lower left of the plot are somewhat unrealistic since low spatial correlation (h 0 ) and a low variance on the prior ( 2 0 ) will not provide a reasonable enough range around the prior to provide a realistic posterior sink estimate which generally is thought to range between 0 and 1.5 PgC/yr (Schimel et al., 2000;Gurney et al., 2002) inter-annually. Increasing either the variance multiplier (along x-axis) or the spatial decorrelation length scale (along y-axis), or both jointly, increases the error variance around the a priori mean allowing more realistic domain-wide summed posterior flux estimates. Therefore if one "de-weights" these sink estimates occurring in the lower left hand portions of the panels in Fig. 9, the RMSE-weighted sink will likely increase to more than 0.57 PgC/yr.

593
The weight of the grand prior (w) has two effects. First, it constrains solutions 594 back towards the prior, essentially anchoring the Kalman filter so that, over time, it does 595 not drift too far from the prior. Given the fact that this grand prior is fixed in time, it also 596 provides a degree of variance inflation (over the regular KF) by providing a lower bound 597 on the prior variance for each filtering step. It is interesting to note that, for cases in 598 which the global prior is weaker (bottom two panels), the maximum sink estimate occurs 599 on the inside of the plot bounds and not at the boundary. The Kalman filter becomes 600 more entrenched without the grand prior since there is no lower limit on the prior 601 variability at each inversion filter step and there is no inflation. Therefore it is likely that 602 the initial reduction in respiration and associated "sink" of carbon in the early months of 603 the year becomes entrenched and leaves a strong sink signature on the rest of the year 604 resulting in the largest sink estimates. We did not test any additional forms of variance 605  Fig. 9. Sensitivity of (a) sink estimate and (b) root mean squared error to varying covariance parameters in inversion. For example, a set of parameter values like w=0, h 0 =600 km, 0 =0.35 provides the estimate w/the lowest RMSE (b) and an estimated sink of approximately 0.55 PgC/yr (a). Nevertheless, qualitatively, this represents a maximum departure from the prior and thus must be viewed with some skepticism due to the likelihood of overfitting the data.
The weight of the grand prior (w) has two effects. First, it constrains solutions back towards the prior, essentially anchoring the Kalman filter so that, over time, it does not drift too far from the prior. Given the fact that this grand prior is fixed in time, it also provides a degree of variance inflation (over the regular KF) by providing a lower bound on the prior variance for each filtering step. It is interesting to note that, for cases in which the global prior is weaker (bottom two panels), the maximum sink estimate occurs on the inside of the plot bounds and not at the boundary. The Kalman filter becomes more entrenched without the grand prior since there is no lower limit on the prior variability at each inversion filter step and there is no inflation. Therefore it is likely that the initial reduction in respiration and associated "sink" of carbon in the early months of the year becomes entrenched and leaves a strong sink signature on the rest of the year resulting in the largest sink estimates. We did not test any additional forms of variance inflation on the model and acknowledge that additional efforts are needed to construct more robust filter techniques.

Comparison to CarbonTracker flux estimates
Given the fact that the majority of the underlying observations supporting the inversion were also used in the Carbon-Tracker project, one would expect posterior flux estimates to be somewhat similar. One of the most important differences between these inversions and CarbonTracker is the optimization of encompassing global fluxes, which affect CO 2 concentrations within our domain. However, this can be mitigated somewhat by the use of optimized CO 2 concentrations from CarbonTracker in the inversion. Under this scenario, one would expect the inversion results to be similar to CarbonTracker but there are still many differences. As can be seen in Fig. 10, the carbon fluxes in the priors, CASA and SiB3, play an important role in the posterior estimates. The posterior estimates of both inversion models display the signature of the a priori fluxes prominently. These results would lead one to believe that either the data does not provide sufficient constraint or the covariance structure is specified too tightly around the prior.
The results can be aggregated up to SiB biomes ( Fig. 1) which are presented in Table 2. The aggregated results show somewhat close agreement with Peters et al. (2007). Because of the differences in land cover classifications between Peters et al. (2007) and the State of the Carbon Cycle Report (SOCCR, 2007), it is difficult to directly compare results. Furthermore, as was shown earlier, the magnitude of the domain-wide annual NEE sink is very sensitive to inflow assumptions although the spatial configuration of the sources and sinks seem much more robust. In this fashion, the proportion of the annual sink due to forests, both conifers and deciduous, ranges between 30% and 37% depending upon inflow choice. This is certainly within the confidence bounds of the estimate given in the SOCCR report as well as very similar to the estimates given in Peters et al. (2007).
One of the most significant contrasts is the placement of the carbon sink on agricultural lands. Peters et al. (2007) shows a large effective sink in the northern Great Plains centered near the state of Iowa where there are very large expanses of corn fields. Our results indicate a strong effective sink in the southern portion of the Great Plains more in the vicinity of large wheat growing operations. Ongoing research in this region of the United States (http:// www.nacarbon.org/nacp/mci.html) seems to validate the existence of an effective sink in the northern Great Plains while  Top panels concern CarbonTracker and lower panels concern our inversion. Left panels show a priori NEE, middle panels show inversion adjustment, and right panels show a posteriori NEE.
agricultural statistics and other recent work (Riley et al., 2009) seem to validate a possible southern Great Plains effective sink as well.

Comparison to filled level 4 Ameriflux data at Southern Great Plains
Posterior respiration and GPP estimates from the model can also be compared to Ameriflux level 4 data. As indicated earlier, there is a spatial representation mismatch in doing so due to the fact that the model estimate is an average over approximately 1600 km 2 and the associated flux tower estimate is over a much smaller footprint, likely less than 1 km 2 . Nevertheless, some useful comparisons and observations can be made. Figure 11 shows comparisons of the model to the observations for weekly ER and GPP at three Ameriflux sites, which appear in the more observation constrained portion of the model domain. The ARM site is one of the more constrained sites in the domain and lies in a relatively homogenous landscape making it an excellent candidate for analysis.
The prior site NEE estimate appears to be improved on average by the posterior flux estimates. In particular, the prior model is corrected significantly in the summer when it predicts significant respiration occurring. Clearly one can see an early spring winter wheat signal in the observations, forming a significant amount of carbon drawdown over an 8-10 week period. SiB3 necessarily balances GPP and ER annually and is thus forced to redistribute this carbon into respiration in other portions of the year. This is the likely reason for displacement of the prior estimate in the summer. The posterior corrects for a large portion of this but the large distance between the prior and observed fluxes make a complete correction difficult. Just as important, but perhaps more subtle, is the fact that the inversion is able to provide significant corrections to ER and GPP separately. SiB3 appears to significantly overestimate GPP. However, due to the annual NEE balance constraint, SiB3 will overestimate ER as well, providing an NEE signal that appears very reasonable. If the forward model is only compared to NEE estimates at various sites then this fact can be easily overlooked but is likely very important to biosphere dynamics on certain time scales.

Evaluation of annual NEE source/sinks against ancillary data and hypotheses
Using two sets of boundary conditions, we arrived at a final sink estimate of 0.5-0. On the other hand, perhaps the globally based sink estimates are too high. The recently completed SOCCR report provides an inventory-based sink estimate for North America of approximately 0.66 PgC per year (land sink) using a variety of data sources collected over the last ten to fifteen years. Uncertainty is presented as a 95% confidence interval, 0 to 1.32 PgC. This is similar to what we have recovered in these inversions. However, this is a mean sink estimate over many years and 2004 is believed to be a year in which the sink in North America was very strong, likely putting the SOCCR estimate closer to 0.8-0.9 PgC/yr, the upper range of their annual estimates. Stephens et al. (2007) called into question the magnitude of the Northern Hemispheric (and North American) global annual NEE sink which has been a cornerstone of inversion results for the last 10 years (Fan et al., 1998;Gurney et al., 2002;Peters et al., 2007) indicating that it may be much smaller than previously assumed. In any case, the rapid expansion of the calibrated CO 2 tower network (currently over 30 towers in North America) should soon provide significant additional data constraints to researchers performing atmospheric CO 2 based inversions.
The spatial character of the annual NEE estimate has several distinctive features. The most definitive feature of the annual NEE estimate shown in Fig. 7 is the large sink located over Texas, Louisiana, Arkansas, and portions of Oklahoma. This sink is located largely between, and to the east of, the ARM and WKWT sites in south central portion of the domain. At first glance this may appear to be an artifact of incorrect transport, poor boundary conditions, or incorrect fossil fuel emissions specifications. However, summing the ARM NEE observations for the year provides a sink estimate of approximately 275 g/m 2 , similar to the estimates the inversion produces to the south of the ARM site (Fig. 7). A likely hypothesis for this sink is the lateral export of crops, primarily winter wheat that draws most of its carbon from the atmosphere in the spring and then is harvested and exported in early summer. The WKWT tower concentrations have proven to be somewhat difficult to model given a number of factors. CO 2 observations at the top of the tower did not appear to be well mixed until well after 12 p.m. LT.
Additionally, the tower is located relatively closely to both the model boundary and the ocean and is in close proximity to fossil fuel sources of major metropolitan areas and oil refining facilities near Houston and Galveston. The aforementioned sink also extends to the east and northeast of the ARM tower. This is an area of significant crop production, with corn and soybeans being grown extensively in the northern portions while soybeans, rice, and other crops are grown to the south in the Arkansas/Mississippi region. This area is also covered by heavily managed forest regions, which produce large annual harvests of wood primarily for paper pulp. These managed forests are largely composed of very young productive loblolly pine trees providing a major source of carbon sequestration. This area is known for quite variable precipitation patterns and it would seem to reasonable to assume that young productive forests in this area would be very productive under the unusually wet and cool conditions of 2004.
It is interesting to note that the most intensely cultivated portion of the Midwestern United States, centered on the state of Iowa, shows little to no sink. This is an area typically planted extensively with corn, which has been shown to be an extremely effective consumer of atmospheric CO 2 . The a priori estimate of NEE based upon SiB3 included a very strong summer time sink of carbon over the Iowa region using a C 4 photosynthesis scheme from Collatz et al. (1992). Lokupitiya et al. (2009) illustrated that a phenology-based model for the fluxes agricultural crops provided flux estimates much closer to those of eddy covariance towers in the region. Additionally, the area of crop production in the a priori model was very crude and did not match the spatial extent of crops, or the mix of different crop types, as given by United States Department of Agricultural maps. The general consensus, as mentioned by Peters et al. (2007), is that this sink occurs over a relatively small intensively farmed area of the country while the agricultural products produced (the effective sink) are distributed out uniformly over the country, effectively spreading out the associated respiration signal. Emerging research from the Mid-Continent Intensive portion of the North American Carbon Program appears to support this hypothesis to some extent (Tris West, personal communication, 2009). Therefore, it is likely that the "missing" sink in this area can be attributed to a poor agricultural crop prior in SiB3 and/or a lack of CO 2 measurements in the vicinity. Although the existence of a sink is widely supported by agricultural crop statistics, the strength of the sink is currently unknown and its estimation is complicated by a number of factors.
First, annual NEE estimates from the corn-planted Bondville, IL Ameriflux site indicate a sink on the order of 500-600 g/m 2 . Soybeans can be expected to provide sinks of about half of this. Assuming steady state conditions over several years, these types of sinks can be attributed directly to the harvest. Approximately 20% of the corn harvest and 35% of the soy harvest is exported overseas, mostly for animal feed, while half of the corn and soy retained in the United States is used to feed livestock domestically (National Corn Growers Association website: http: //www.ncga.com/files/pdf/2009WOC.pdf, Soy Stats, http:// www.soystats.com). Most of the carbon in this livestock feed is then returned to the atmosphere as CO 2 and CH 4 at locations where it is consumed by livestock. Almost 70% of the feedlots in the United States are located in just 3 states: Texas, Kansas, and Nebraska (http://www. cattlenetwork.com). This may provide a partial explanation for the lack of an agriculturally-induced sink over Nebraska and Kansas, states with very high crop production and intense livestock operations, and the existence of sinks over portions of Arkansas, Mississippi, Missouri and Illinois, states with relatively high crop production but with significantly less livestock operations. A rapidly evolving ethanol industry in the area further complicates the picture. Currently, this is an area of intense research (http: //www.nacarbon.org/nacp/mci.html)and one may expect a much more complete picture to emerge concerning the carbon balance of the Midwest United States within a few years. An important point to keep in mind is that the addition of a carbon sink in the Midwest United States would likely be correlated with weakened sinks (or increased respiration) in other areas of the domain in order to constrain the annual domain wide source/sink estimate.
Forested regions in the northwestern United States and boreal forests of Canada show slight sinks. However, variability estimates surrounding these sink estimates are typically much smaller than the variability estimates of similar sink magnitudes in the Midwest or southeastern United States showing relatively more confidence in the sink despite the lack of proximity to the observing towers. The sink estimate in the northwestern United States is not surprising since the northwestern coastal mountains of California, Oregon and Washington have been intensely managed over the last 50 years and are expected to provide a sink of carbon for many decades into the future (Alig et al., 2006). The estimate for the boreal forest regions appears much harder to objectively evaluate. Most studies have indicated that Canadian ecosystems should currently be a weak sink, although the projection of this weak sink into the future is highly uncertain. The inversion results show a fairly carbon neutral Canada on average, but shows the boreal forests of central Canada and the boreal and coastal forests of western Canada as slight sinks while the agricultural plains of Canada and the forests of eastern Canada provide slight sources. It is interesting to note that areas to the south of the two Canadian towers show an annual source of carbon in an area just to the east of large expansive forest ecosystems of British Columbia that have recently experienced unprecedented bark beetle invasions and tree mortality. It is important to note that forest fires were not included in the SiB3 domain run for the regional inversion. Average carbon emissions from Canadian forest fires were estimated at 27±6 Tg/yr (Amiro et al., 2001), a non-trivial amount that could increase the strength of the boreal forest sink predicted by the inversion.

Conclusions
GPP, ER, and NEE flux corrections implied by this inversion provide posterior annual NEE estimates similar to those provided by a number of independently derived models including CASA (via CarbonTracker optimized) and the MODIS 17 GPP product. NEE estimates for the entire domain appear on the low side of estimates derived from global models, which is understandable given the lack of constraint on some key regions of high annual GPP, and hence potentially high annual NEE. This was corroborated by a comparison to INTEX aircraft data which shows the existence of a deficit in GPP over the southeast which would, when all other things are considered equal, inflate the domain-wide sink closer to levels estimated from global models such as CarbonTracker. Results are relatively sensitive to a number of parameters in the inversion setup, which is also to be expected with an inversion constrained by such a sparse observing network. Using a temporally uniform boundary condition seems to produce a very unrealistic annual sink on the order of 0.38 Pg per year, supporting the notion that regional inversions require realistic boundary inflow of CO 2 . However, much to our surprise, we find that two completely independent boundary inflow estimates provide similar results with the main difference being an approximately 30% difference in magnitude. This leads us to believe that, while probably not preferable to optimized global CO 2 fields, the inclusion of annual NEE balanced models (such as SiB3) in global models used to provide boundary inflow estimation does not significantly damage inversions based upon it.
In the course of trying to improve NEE estimates, we were able to find that the inversion was able to provide some degree of correction to the individual summands of NEE, ER and GPP, which are generally highly correlated at many different scales in time and space. Considering that SiB3 currently calculates ER as a relatively simple function of soil moisture and temperature such that annual ER equals annual GPP, the significant adjustment inferred upon GPP may prove to be valuable estimation of other quantities of interest in the biosphere. For example, while photosynthesizing, plants must generally release water to compensate, meaning that artificially high GPP may infer artificially high water exchange with the atmosphere and possibly associated latent heat fluxes.
The agricultural Midwestern United States appears to play a large role in the inversion results, providing a large sink. However, the sink does not correlate exactly with crop productivity, when compared to crop production maps from the United States Department of Agriculture, and several states with significant crop production such as Nebraska, Kansas, and Iowa, appear to be in approximate annual carbon balance. While the magnitude of this difference between carbon neutral states with crops and carbon sink states with crops is likely influenced by the lack of data in the inversion and the general unconstrained nature of the solution at fine scales, the discrimination between them seems likely to stay. One hypothesis proposed is the lateral movement of crops which has been shown to be a major portion of the carbon budget globally . The main crops of interest in the domain are wheat, soy and corn. Soy and corn are grown across large expanses of the north-central Midwest and are primarily used to feed livestock. These livestock are typically fed in feedlots in the states of Iowa, Colorado, Nebraska, Kansas, and Texas, generally located to the west and south of the areas of growth and harvest. The end result would be that eastern states within the Midwest would be a sink because of the near complete export of crops grown there. However, states in the western portion of the Midwest would receive the majority of these crops where they would be fed to cattle and other animals, returned to the atmosphere as CO 2 and CH 4 and largely balance any local sinks due to crop production.
Technical considerations concerning the inversion could also affect these results. In particular, a large amount of missing data for the WKWT (Moody, TX) tower leaves the southern boundary inflow unconstrained beyond the normal PCTM inflow. This could result in the inflation of an Oklahoma/Texas sink to account for a positive bias in the inflow at the southern boundary, particularly after 1 July 2004 when the Midwest receives its heaviest influence from the Gulf of Mexico. The WLEF tower was also missing most of its observations for June, a time of intense drawdown for croplands to the south of the site.
In 2004, the southern states of Texas, Oklahoma, Kansas, Louisiana, Arkansas, and Mississippi had an extremely wet summer, potentially mitigating some degree of drought and providing an increase in GPP for the region which includes managed forests, a large percentage of the United States' exported wheat crop, and soybeans and other crops along the lower Mississippi river valley. Additional research is needed to determine if any of these could represent a plausible hypothesis that would result in the net carbon neutrality of large crop growing states in the western portions of the Great Plains and the expansive southern and Mississippi river valley sink predicted by the inversion.