I think the authors have developed a very thoughtful inverse modeling system based on hierarchical Bayesian statistics. The paper is well-written, the statistics are carefully explained, and the model is thoroughly evaluated in the manuscript. I think the inclusion of observational biases, correlated errors, and the calculation of posterior uncertainties are major strengths of the proposed modeling framework. I highly recommend this article for publication.

My one concern is the use of relatively coarse basis functions within the inverse model. On one hand, I understand the necessity of reducing the number of unknowns in the inverse problem. Without the dimension reduction, I doubt it would be computationally possible to estimate the fluxes using an MCMC algorithm (not to mention the number of CTM runs that would be required). On the other hand, I suspect that the quality of the flux estimate will depend on the number and quality of basis functions used in the inverse model. The basis functions used in this study are for relatively large TransCom regions and vary by month, so the inverse modeling framework proposed here would not be able to resolve patterns in fluxes at smaller spatiotemporal scales. This may not be an enormous concern for simulations using OCO-2 but could be a real challenge if using observations that are more sensitive to finer-scale variability in CO₂ fluxes (e.g., aircraft observations or observations from future satellites). In addition, existing bottom-up models of biospheric CO₂ fluxes yield extremely different spatial patterns (e.g., see the MsTMIP model inter-comparison project), and I suspect that the flux estimate might look different if the basis functions were drawn from a biospheric flux model with very different spatiotemporal patterns from CASA. At a more technical level, the basis functions have different spatial support, and I wonder if that could/should impact the hyperparameters on alpha (particularly when the hyperparameters are informative). Also, the spatiotemporal patterns of the fluxes within each basis function are assumed known, and I wonder if that assumption could reduce the size of the posterior uncertainty bounds in a way that is undesirable. Overall, I think the inverse model is very thoughtfully-constructed, and I don't want to suggest any major changes. With that said, I think the authors may want to head off or preempt some of these concerns within the manuscript. While it is true that some inverse modelers use a 'scaling factor' approach, other inverse modelers do not for the reasons outlined above. Ways to head off these concerns could include (1) provide additional guidance on how to choose basis functions, particularly in light of uncertainties in existing biospheric CO2 flux estimates; (2) discuss the computational feasibility of using a greater number of basis functions (e.g., with an eye toward future studies), and/or (3) include more discussion on the trade-offs between using a greater number of basis functions versus the reduced computational demands of using fewer basis functions. Overall, I think anything that can be done to head off these concerns within the manuscript would strengthen the manuscript.

Specific suggestions:

• Introduction: For many statements in the introduction, there are numerous different citations that you could cite to support the claim. In these cases, you might want to preface the parenthetical citation with "e.g.,". That way, you don't imply that the cited reference is the only reference available that supports this claim. Rather, it is merely one reference, among other possible references.
• Line 53: Most Bayesian-synthesis inversions use prior distributions to encode uncertainty in the fluxes. Can you be more specific about the innovation in WOMBAT? E.g., do you add in a hyperprior, or some other novel feature?
• Lines 58-60: There are a handful of existing papers that have used Monte Carlo simulations or conditional realizations to estimate posterior uncertainties in global trace gas inversions (e.g., papers by Frederic Chevallier and Junjie Liu). However, these papers (to my knowledge) do not use a full-blown MCMC framework nor do they propagate uncertainties in the hyperparameters.
• Line 82: Michalak et al. (2004) did not set the prior flux field to zero everywhere. Rather, the prior flux field was non-informative with an unknown mean that is solved as part of the inverse model.
• Line 101: Should there be a comma after "WOMBAT"?
• Sect. 2.2: It feels slightly confusing to use the variable "Y" to denote both fluxes and atmospheric mole fractions.
• Line 172: I believe the column average is also based on the pressure weighting function (not just the averaging kernel).
• Line 202: What do you mean by "statistical efficiency"?
• Line 225: Michalak et al. also estimated uncertainties in the hyperparameters but did not explicitly integrate those uncertainties in an inverse model. They also did not use a hyperprior. In my experience, some parameters are relatively easy to estimate using maximum-likelihood (ML) (e.g., the variances of the covariance matrices) while others are challenging to estimate from the atmospheric observations using ML (e.g., the covariance structure). I.e., often there isn't sufficient information in the atmospheric observations to guide the choice of these parameters. I think one advantage of using a hyperprior is that there is prior information to inform the hyperparameters even when the atmospheric observations are not sufficient to inform those parameters.
• Lines 369-370: I don't necessarily agree with the statement that the space-time patterns in existing bottom-up estimates are realistic. Bottom-up biospheric flux estimates from recent inter-comparison projects show widely differing spatial and seasonal patterns. For example, there are large spatial and seasonal differences among models participating in both MsTMIP and the TRENDY model comparisons. The large differences among existing bottom-up biospheric flux models implies that regional space-time CO₂ flux patterns are not well known.
• Line 370: I wouldn't say that this assumption is ubiquitous. While some inverse modeling systems do make this assumption, there are also plenty of inverse modeling systems that do not make this assumption. For example, inverse modeling systems like the CMS-Flux system at NASA JPL and CarbonTracker-Lagrange system at the NOAA Global Monitoring Laboratory do not make this assumption.
• Line 465: The word "embarrassingly" feels a bit too colloquial here.
• Fig. 3 caption (and other figures throughout the manuscript): I think it would be useful to provide the main takeaway message of the figure within the caption. I think that would help guide the reader on how to interpret the results.
• Line 558: I think there's a typo at the beginning of this line.
• Line 582: I would use a more specific term than "deeper" to avoid confusion here. (This wording is also used later in the manuscript.)
• Fig. 4: The shaded areas on the right hand side of this plot are difficult to distinguish from one another. I wonder if there is a way to improve this aspect of the visualization. This suggestion applies to a few other similar figures as well.
• Overall: I think it would be nice to have an appendix with all of the variable definitions, so the reader can keep track of what's what without searching through the manuscript text.

This paper is very well presented and interesting and contributes significantly to the inversion community.

I support publication after the following concerns are addressed:

1) Supplement vs appendix

The appendix are too abundant. Please keep only Appendix A-B and put the rest in supplement

2) Better identify advantage of WOMBAT

Even though the OSSE and MIP comparison are interesting, it is hard fully separate the impact of WOMBAT compared to classical inversions and part of the patterns are necessarily linked to the CTM used. Using another CTM is beyond the scope of the manuscript, but I would suggest showing an analytical inversion with the same base functions, but no estimation of hyperparameters.

It would show how a classical inversion would perform with GEOS-Chem compared to WOMBAT, which would partly disentangle different error sources

3) Information content

The author should further elaborate on the information content in the system. In particular, contrary to classical inversions, WOMBAT must estimate hyperparameters, which are much more numerous than in standard inversions.

Thus, is there enough observations? When adding new regions is conflicting with the number of observations? What is the cost of one MCMC with the present configuration, and when does it get too expensive to compute with increasing dimension?