Summary and assessment:

This study applies "functional ANOVA", a version of the analysis of variance (ANOVA) method that estimates spatially-varying effects, to CO2 flux estimates from several inversion studies. The work builds upon previous studies, and uses some computationally efficient approximations to covariance structures in order to increase the number of contrasts that can be estimated. It is more of a methodological study, and illustrates how some of the parameters that can be estimated from this methodology can provide practical insights into the findings (i.e., messages that may be relevant for a broader audience than the scientific community). It does not offer major methodological development, but more of a case of demonstrating how this method, with the computationally efficient spatial covariance model, can be used for this application. The writing is of a very high standard, with very few corrections necessary. It is of an appropriate standard for GMD and the manuscript sits well within the Journal's scope.

I believe that this is a well-conducted piece of work, and I would like to see it published. However, there were a small number of issues that I believe need to be addressed, either by some additional analysis or further discussion in the text. I am recommending this as a Major Revision, although it is really on the border between Major and Minor Revisions.

Major comments:

 - I'm concerned that one of the assumptions about the error distribution, namely that the realisations are independent, may not be valid. The realisations are taken as successive months within a three-month block. In my experience, geophysical fields and their associated error distributions are typically correlated in space and time. The authors need to address this issue directly, perhaps demonstrating numerically that the residual fields cannot be distinguished from random draws, or arguing persuasively why this should be overlooked for the sake of illustrative example. One thing to consider here, and probably should be raised in the revision, was whether the inversions were conducted as a "long" run (i.e., at least spanning all of the three months considered), or whether they were run separately for the individual months in question.

 - It is unclear what is going on with the ocean areas. The distribution of fluxes is very different for land areas or oceans. Similarly, the distribution of fluxes depends very much on the type of vegetation cover. Consider for a moment differences between tropical rainforest vs savannah vs deserts vs high-latitude boreal forests. Are we assuming a uniform standard deviation? This sould be discussed. 

 - Also on the topic of oceans, it seems that the fluxes for most of the oceans are zero, or close to zero on the colour scales shown, except for some small areas (e.g., over the Mediterranean or Red Seas in Figure 10). Is this intentional? Is a land-mask being imposed? I recommend discussing this further.

Minor comments:

 - Multiple places: There are multiple references to the supplement of this manuscript. I recommend making the references more specific, such as to particular figures, tables, sections, pages, etc.

 - Abstract: I thought this was rather short. I recommend adding a few sentences expanding upon conclusions and implications of this work.

 - L3: "are combined" -> "can be combined"

 - L4: "impact flux estimates" -> "impact such flux estimates"

 - L8: I found the term "mode differences" somewhat cryptic. Please review and reformulate.

 - L31: The term "data fusion" isn't particularly well defined. I suggest adding some clarifying remarks or using a different term.

 - L37: "ANOVA methodology" -> "The ANOVA methodology"

 - L41: "multiple component" -> "multi-component"

 - L47-48: Is the method of Cressie et al. (2022) based on functional ANOVA? If so, I suggest noting this.

 - Page 2: Reading this far, I was unsure whether the authors were talking about gridded data or not. This became clear later on, but I suggest making a note about this somewhere around this point in the Introduction.

 - Page 2: I was left wondering how functional ANOVA handles different spatiotemporal resolutions (or, more generally, different spatiotemporal partitioning). Do all of the data need to be on a common grid? 

 - L63-64: References to Sections 2 and 3 are swapped around. Review.

 - L69: I suggest noting the units of XCO$_2$

 - L70: Do you ave a version number for the ACOS retrieval algorithm?

 - L73: "within a single polar orbit" - do you mean "within a single model grid-cell"?

 - L75: "a precise" -> "a more precise"

 - L83: I find URLs within the text, particularly when it is messy such as DOIs, rather unsightly. I suggest creating a bibliography entry for this dataset, and shifting the URL there. 

 - L101: "SSDF" - Has this acronym been introduced? 

 - L103: "The variable of interest" -> "The variables of interest"

 - L122: "the individual retrievals are both variable and moderately correlated in space and time" -> "the individual retrievals are both uncertain and the associated errors are moderately correlated in space and time"

 - L122: What is meant by "moderately" here?

 - Figure 1 caption: I don't think the acronym "JJA" has been introduced until now.

 - L123: "spatially aggregated" - as a gridded or regional average?

 - L129: "have similar resolution" -> "have similar spatial resolution"

 - Figure 2 caption: "over North America." -> "over North America, with columns for the different inversions."

 - Table 1: I recommend adding an extra column with a key reference for each retrieval or each retrieval dataset.

 - L166-167: Here we see the issue of months within a season being used for different replicates (see the first of my "Major comments" above).

 - L170: "are often" - why not "must be"? Is there an alternative?

 - L172: "the factor and interaction effects add to zero, e.g." -> "the factor effects add to zero, i.e."

 - Equations below L172, and L175: Why not put these on the same line? Would they not fit?

 - L195: Here the realizations are assumed to be independent (see the first of my "Major comments" above).

 - L214: $\mathbf{\theta}$ was not defined, as far as I could see.

 - L215: What do you mean by "diffuse"?

 - L223: "for the component and the noise" -> "for the components and the noise"

 - Page 11: Here I was hoping to see more about the question of a land mask, or different distributions (e.g., heteroscedastic) for land/ocean areas.

 - Page 12: It gets mentioned much later in the Results section (L285) that negative anomalies correspond to carbon uptake by the land, and positive anomalies to the release of carbon into the atmosphere. I think this is an important part of the interpretation of the results shown, and as such suggest introducing this much earlier in the Results section (e.g., within the first few paragraphs, maybe also in the caption for Figure 4).

 - L240-245: I would have like to see a bit more interpretation. For example, what are the implications of the estimates for the range parameter being an order of magnitude larger for the mean and year effect than the other terms? What does it mean that the error standard deviation is so much larger than that of the other terms? What can one read from the relatively narrower credible intervals for the error standard deviation than those of the other terms? 

 - L254: "The left panel shows Pr(|µ(s)| > |β*(s)|)" - How was that estimated? With reference to the posterior samples?

 - L260: It would be worth noting in the text that the posterior probabilities for the secondary question, Pr(|α(s)| > |β*(s)|), are clearly lower than for the primary question, being Pr(|µ(s)| > |β*(s)|).

 - L285: "negative anomalies (increased uptake)" - as noted above, this should be mentioned earlier.

 - Figures 7 & 10, bottom panels (data source effect): I would have been interested to see the impact of the locations of the in-situ sites. Consider adding dots for the locations of the CO2 monitoring stations.

 - L295: "The prior fluxes over the continent (not shown)" - Why not put this in the supplement?

 - Around L293-294: I would suggest noting that for much of the domain, the direction of the flux estimated by the posterior is ambiguous.

 - Last paragraph on Page 15: I would recommend expanding some of the interpretation. On the question of the data sources, adding to the Figure 10 the locations of the monitoring sites could provide some basis for some of the interpretation of the data source effect.

 - L320: "The functional ANOVA identified local consensus flux anomalies for both continents in the presence of variability across inversion systems and atmospheric CO2 data sources. " -> "The functional ANOVA identified local consensus in flux anomalies for both continents across different inversion systems and atmospheric CO2 data sources. "

 - L339: "heterogeneous across space" - I think there's a word missing here. Do you mean heterogeneous standard deviations?

 - Discussion: I think it would be useful to expand upon the importance of spatial heterogeneity. See the second of my Major comments. 

 - References, multiple places: CO2 should be CO$_2$ and XCO2 should be XCO$_2$

 - References, multiple places: at least two of the author lists are truncated (finishing with "et al.") while others extend to many authors. Please make this consistent, preferably checking the GMD Guide to Authors.