Interactive comment on “Geostatistical inverse modeling with very large datasets: an example from the OCO-2 satellite”

We would ﬁrst like to thank Ian Enting for his thoughtful review. His input and sugges-tions have been extremely helpful, and the manuscript is better for it. Below, we have listed Dr. Enting’s comments in bold and discuss the associated changes we have made to the manuscript. We would like to thank Peter Rayner for his feedback and constructive ideas on the manuscript. These ideas have been very helpful for improving the overall quality of the manuscript and strength of the analysis. We have listed Dr. Rayner’s comments in bold typeface below and discuss the associated changes we have made to the manuscript. Abstract. Geostatistical inverse modeling (GIM) has become a common approach to estimating greenhouse gas ﬂuxes at the Earth’s surface using atmospheric observations. GIMs are unique relative to other commonly-used approaches because they do not require a single emissions inventory or a bottom-up model to serve as an initial guess of the ﬂuxes. Instead, a modeler can incorporate a wide range of environmental, economic, and/or land use data to estimate the ﬂuxes. Traditionally, GIMs have been paired with in situ observations that number in the thousands or tens of thousands. However, the number of available 5 atmospheric greenhouse gas observations has been increasing enormously as the number of satellites, airborne measurement campaigns, and in situ monitoring stations continues to increase. This era of proliﬁc greenhouse gas observations presents computational and statistical challenges for inverse modeling frameworks that have traditionally been paired with a limited number of in situ monitoring sites. In this article, we discuss the challenges of estimating greenhouse gas ﬂuxes using large atmospheric datasets with a particular focus on GIMs. We subsequently discuss several strategies for estimating the ﬂuxes and 10 quantifying uncertainties, strategies that are adapted from hydrology, applied math, or other academic ﬁelds and are compatible with a wide variety of atmospheric models. We further evaluate the accuracy and computational burden of each strategy using : a ::::::: synthetic : CO 2 observations from :::: case ::::: study :::: based ::::: upon NASA’s Orbiting Carbon Observatory 2 (OCO-2) satellite. Speciﬁcally, we simultaneously estimate a full year of 3-hourly CO 2 ﬂuxes across North America in one case study – a total of 9 . 4 ⇥ 10 6 unknown ﬂuxes using 9 . 9 ⇥ 10 4 ::::::: synthetic observations. The strategies discussed here provide accurate estimates of CO 2 ﬂuxes 15 that are comparable to ﬂuxes calculated directly or analytically. We are also able to approximate posterior uncertainties in the ﬂuxes, but these approximation are typically an over- or underestimate depending upon the strategy employed and the degree of approximation required to make the calculations manageable.

observations. We have modified the text in the abstract to clearly indicate that these are synthetic observations.
• A complicating aspect is that the fluxes from carbon-tracker are themselves the product of an inverse calculation and so will have different spatiotemporal correlations than actual fluxes.
This is also a fair point. We do not know what actual CO 2 fluxes are, but a flux product like CarbonTracker is likely the closest approximation to actual fluxes that one can find. To this end, we have added a caveat on pg. 11, line 25 (of the original GMDD manuscript); real-world CO 2 fluxes may be different from the CarbonTracker fluxes used here. Hence, the analysis presented in this study is an approximation or prototypical example of the computational challenges one might encounter in an inverse model.
• I would strongly disagree with the claims of uniqueness of GIM with regard to including other types of information.
There are two different things: 1. The geostatistical approach of using spatio-temporal correlation structure as a technique for regularising an ill-conditioned inverse problem; and

The inclusion of additional information about fluxes (incorporated into GIM through β ).
There is nothing to prevent the inclusion of additional information into inversion techniques that do not use geostatistical constraints.
While the specific form of p(s|β) that is used leads to linear equations and a direct solution, once direct solutions are replaced by 'variational' approaches, more general forms of p(s|β) can be incorporated, either with or C2 without the use of regularisation by imposing a spatio-temporal correlation structure.
We agree with point #1. Recently, numerous inverse modeling studies that use a classical Bayesian approach have included spatio-temporal correlation structure in the prior covariance matrix (referred to in this manuscript as Q). Hence, while some previous GIM studies have argued this is a unique feature of a GIM, we do not make that argument in the present manuscript. We have modified line 1 on page 4 to make this distinction clearer.
With regard to point #2: we agree that more general forms of p(s|β) could be incorporated, but this has not been done in existing atmospheric inverse modeling studies. I.e., there is no atmospheric inverse modeling study that we are aware of that has done that to date. In theory, one could design an inverse model in any number of ways with different distributional assumptions, prior probability densities, hyperparameters, and hyperpriors. However, nearly all atmospheric inverse modeling studies to date use a relatively narrow set of cost functions that have a similar form. In that context, the way in which GIMs have been applied to atmospheric inverse problems and the way in which these studies have formulated the prior probability is unique relative to what has been done to date in classical Bayesian inverse modeling studies using atmospheric trace gas observations.
• page 11, Line 5. The Miller et al. (2018) study doesn't seem to provide much information about the actual spatio-temporal correlation structure of the OCO-2 data (i.e. the structure of R. More discussion of this would be desirable.
We have added several sentences to the revised manuscript to clarify the inverse modeling approach used in the manuscript.
The actual error correlation length and correlation time varies depending upon the region in question and the day of the year, depending upon factors like ob-C3 servation type (e.g., nadir versus glint), atmospheric aerosol concentrations, and variations in surface albedo (e.g., O'Dell et al. 2018). In Miller et al. (2018Miller et al. ( , 2019, we estimated a spatial error correlation lengthscale between 0km and 3,800km for land nadir and land glint observations (version 9) using a spherical covariance model. We also calculated error correlation times between 0 days and 45 days, also using a spherical model.
In the present manuscript, we used a diagonal structure for R. All existing inverse modeling studies that we are aware of using OCO-2 observations utilize a diagonal structure for R, and we use a similar approach that is prototypcal of existing inverse modeling studies. Also note that the STILT footprints here were run every 2 seconds (approximately every 15km) along the OCO-2 flight track, so the synthetic observations used in this study are spaced further apart than the observations in the full OCO-2 dataset. In some locations, this spacing is comparable to the estimated spatial error correlation length in the OCO-2 observations.
It may be helpful to incorporate off-diagonal elements within R for future inverse modeling studies -to better account for the information content of the OCO-2 observations and/or to estimate rigorous posterior uncertainties. However, that has not been done in existing OCO-2 studies to date, and doing so would likely necessitate accounting for the highly non-stationary structure of these off-diagonal elements. We discuss the possibility of incorporating off-diagonal elements within R in the manuscript (e.g., pg. 13 of the GMDD manuscript). The overall focus of this paper is on computational approaches to large inverse problems, and to that end, we felt that developing an approach to describe non-stationary error covariances in R was beyond the scope of the current manuscript.