Review of "Predicting and correcting the influence of boundary conditions in

regional inverse analyses", Nesser et al., for GMD

This manuscript is well-prepared, well-organized, and well-written. It is also a very useful study for practitioners of atmospheric inversions at regional and local spatial scales, where background errors are complex and often large in magnitude, because decreasing spatial scale usually increases variability (and error) in the background estimate. Analysis at these scales is becoming more and more common, as researchers seek to provide emissions data at "actionable" scales, i.e. for a given oil and gas basin, city, or state. This should be published with minor edits.

Most of my comments below are requesting clarifications of legends/figures and of some more intuitive explanation of one of the results -- although I believe part of my confusion was due to a mis-labeled figure.

Detailed (minor) comments:

L16 "a gridded flux estimates" should be made either singular or plural

L196 two periods after p

Fig S1: Having trouble understanding the lines in the lower three panels relative to the legend - both the "no correction method" and the "downwind error" are in a similar shade of gray, but I do not see two different lines in the plots.  Also what is the light blue dashed line (e.g., around 0.2 in panel (d))?  I wonder if the color palette could also be changed to show larger differences between the lines, while still making them color-blind friendly?

What about the boundary method, doesn't it reduce upwind differently from downwind?

L317: Related to the above is the understanding of the text here - I am not seeing the boundary method reducing both upwind and downwind error in panel S1(d), I only see one line there.  Also the buffer method does seem to reduce downwind errors relative to the "no correction method" case, I think?  I think I am just not understanding the legend, so please clarify.

Fig. 4 - perhaps title of the left panel could be changed to "information ratio" for consistency.

L352 should say "outperforms the buffer method in all cases"

L355 -- is this likely to be the case for other inversions?  (that the preview underestimates the influence length scale with no correction by ~50% but overestimates it for the boundary adjustment by ~50%).  I imagine that would be informative for people who would like to use it.

L364 - is the tower based non-IMI inversion included in these refs, perhaps it should be.

Refs Varon et al 2023a and 2023b are the same

L375 and Fig 5.  I believe (d) and (e) are swapped when referred to in the text.  This led to some confusion on my part but the overall question is why does the preview give such different spatial results from the real error?  (i.e. why are these two panels so different)? Perhaps some intuitive explanation could be made here, and explanation with what regime we are in in terms of R.

L387  - please clarify - the boundary method reduces computational cost relative to the buffer method, but still increases it relative to no correction ,yes?  How much more expensive is the buffer method in this example? How many elements are being estimated in the boundary method (one value per observation, or one value per edge)?

L395-400 It seems the preview-derived uncertainty spatial map is not correlated (or maybe it is negatively-correlated) with the map of the prior (and prior error). But the prior uncertainty does factor in through the beta term, can this be explained?

L425 Again clarify: reduces computational cost relative to the buffer method.

L445 - this was a useful summary of the difference in the preview vs. the full error, answering some of my questions on the earlier section, but I still think some reader confusion could be alleviated with an intuitive interpretation of the results in the previous section.

Supplement notes:

S1, line 8 add the definiton of xj:  "The area A converts the flux xj from units of mass per area...".

S2: derivation of Eq. 9.  Says we will show that Eq. 10 is equal to delta xhat = -Gec in the extremes of R, this should read Eq. 9 I believe?

Fig S1 as noted earlier - I had trouble interpreting that figure due to the choice of colors and legend.

This work examines the role of boundary condition (BC) errors on regional inversions of long-lived trace gas emissions. As summarized in the conclusions, the findings of this work are preview (prediction) and diagnostic metrics of BC induced errors in the posterior fluxes, and the comparison of two inversion methods for treating BC errors.  These are derived and examined through 1D, two-box, and 2D examples. While the aims of this work are laudable, overall, much of it struggles to demonstrate relevance, accuracy, or practical usefulness to the problem at hand.  In particular, the applicability of the 1D and 2 box results,  and the preview estimates is not clear.  Several of the 1D results don’t translate to the findings of the 2D case, in terms of the dependence of the BC error on length from the boundary for any method, or the accuracy of the preview estimate. The preview estimate seems to miss key functional relationships and lacks applicability as a bounding estimate in the 1d and 2-box model; in the 2D application it simply seems inaccurate. 

The diagnostic of BC errors is also presented as a conclusion.  However, this diagnostic (equation 5) is not a new result (it is a straight forward application of theory from Rodgers).  It is also rather useless in practice.  While it does allow for correct quantification of the BC induced error on the posterior emissions, it requires precise knowledge of the boundary condition error (epsilon_c) to begin with.  In practice, if one does know the BC error to begin with, one corrects the BC’s (or the obs) by this amount before doing the inversion.  More useful diagnostics would be knowing how errors in an estimate of the BC bias may impact the solution, upper bounds on the impact of BC errors under particular conditions, etc.  

Thus, the main useful findings of this work are related to the comparison of the buffer method (clustered or non-clustered) to the boundary method, and I suggest a revised manuscript focus much more on these.  There are some outstanding questions related to their setup and findings on these aspects, see detailed comments below. The application / relevance of the 1D and two-box model results and preview equation need to be further demonstrated prior to their inclusion in a final manuscript. If omitting these aspects, a lot of my detailed comments below can be ignored. 

Comments: 

That the preview estimates generate a “stable estimate of the influence length scale” is an overly generous assessment. In Fig 3, the preview underestimates the length scale over which BC induced error is significant by several a factor of two or more, hitting 0 at 3-4 grid cells while the actual error gets to zero around 7-10 grid cells.  In Fig 4, the preview completely misses the functional dependence on R or tau, and it unfortunately also underestimate the length scale and as such does not provide a useful upper bound.   In the 2D tests of Fig 5 / section 5, the preview of the boundary condition induced error (Fig 5d) looks nothing like the actual boundary condition induced error (Fig 5e). Line 445 claims the preview estimate is accurate near the boundary, but this isn’t rigorously shown, nor is it discussed when / why the preview differs substantially from the actual error. 

Thus, what is the value of the preview? Is it a bound in any way? Does it inform an inversion setup in any way?  The 1D analysis in sections 2 and 3 imply there is a distance from the edge at which point the boundary condition error becomes minimal, usually around 8-10 grid cells.  The preview estimate again estimates this to be the case.  But the actual error shown in e is spread throughout the domain. Thus I’m not sure what is the value of the analysis in sections 2 and 3, which seemed to imply that there was a “safe” distance from the edge even without any boundary or buffer method applied.  The error in 5e also shows a lot of structure connected to the prior, rather than the boundary conditions. I’m thus skeptical about it being “boundary condition induced error.”  The setup also seems to violate a key assumption of the methods, in that x_A doesn’t not seem to be a sample of an emissions distribution with a mean of x_T.  The true emissions x_T are zero in large parts of the domain.  A random sample around this with standard deviation of sigma_A would lead to negative prior emissions.  However, there are no negative prior emissions.  It seems they’ve taken the true inventory and the prior inventory from independent sources.  There’s an inflation to the prior to make some statistics match (Table 1, comment d), but this is rather vague and I’m rather skeptical that it makes x_A a sample of a gaussian distribution with mean x_T and variance s_A.

2D Boundary method: Optimize how many separate elements along the boundary? One for each horizontal and vertical grid cell? In that sense, how is this different than a setup with a buffer method, with a buffer that is one grid cell thick? The only hint I see here is line 104 that says “one or more.”  Given this is the best approach being put forward, some more details could be provided. 

General: I’m a bit confused about what model is being used for which test / set of results.  Let me see if I have this correct. Section 2 presents a 1D (Section 2.2) and then two-box (Section 2.3) model.  Fig 2 shows results from the 1D model.  Fig 3 also presents results from the 1D model, but described in section 3. Why one of these is in the methods and the other the results is a bit confusing.  Also, the text on Fig 3 (line 287) refers to Eq 14, which is for the 2 box model. Also I think they didn’t mean to refer to Eq 14 at all here, since it doesn’t contain epsilon_c — perhaps they meant Eq 13?  Fig 4 caption says it is based on the 1D model, so I would think the “preview” line corresponds to Eq 9.  However the text (lines 343, 347) refers to Eq 9 and Eq 15, and Eq 15 is for the two-box model, not the 1D model.  How is the 2 box model used for the results shown in Fig 4? 

Fig 4: I’m having trouble thinking about how these results are applicable, given that the metric is grid cell number (which seems arbitrary).  Can this instead be presented in terms of a length scale related to wind speed, grid size, and chemical lifetime? The latter would be much more useful and general.  Otherwise, this plot tells me I can get by with a buffer of 4-8 grid cells, regardless of whether the grid cells are 50 m and I’m inverting for CO2 fluxes under windy conditions, or the grid cells are 500 km and I’m inverting for NOx fluxes under stagnant conditions, which doesn’t seem correct. 

427: The claim that this applies to short-lived gases seems to be an overstretch. As alluded to in comments above, I doubt the direct applicability of some of the assumptions and equations developed here for short-lived gases, for whom chemical lifetime needs to be accounted for somewhere, and also given that atmospheric chemistry often leads to nonlinearities in short-lived species that are not accounted for in the present framework.  

26: Does this only apply to long-lived trace gases? It’s not an issue for e.g. NH3, NO2, etc.  Also, regional inversions of trace can concentrations don’t necessarily employ a CTM.  They may use online models, or Lagrangian back trajectory models, plume models,…

70: I don’t agree that a long spin-up makes initial conditions unbiased compared to the BCs.  If the BCs are biased owing to transport, but initial conditions are biased owing to incorrect fluxes within the domain, these biases could be different. 

71: c would include contributions from the ICs as well.  These could be different than those from the BCs, see above.  If one wants to avoid this, one cold assume the inversion period is sufficiently long enough that the contribution from ICs becomes negligible compared to BCs and fluxes within the domain. 

85: This might be a bit clearer if defining epsilon directly as epsilon = y - (K x_T + c)

Eq 6: The last term does not seem correct to me, or there is an inconsistency between the text describing the setup and the equations presented. Perhaps not a big deal as this equation never is actually used, as far as I can tell? But it should at least be correct if presented. My derivation is as follows. 

Retrieval with true boundary conditions and an observing system error eps_T that is owing only to things like obs error, representational error, model error, etc., is (focusing only on the G eps terms for simplicity):

eps_T = y - (K x_T + c_T) 

x_1 = x_A +… + G eps_T

Retrieval with boundary correction method and boundary error eps_c, omitting augmented state vector elements:

eps’ = y - (K x_T + c_T + eps_c)  = eps_T + eps_C

x_2 = x_A +… + G’ eps’ 

= x_A +… + G’ (eps_T + eps_C)

Taking the difference:

x_2 - x_1 = … G’ eps_T + G’ eps_C - G eps_T 

= … (G’ - G) eps_T + G’ eps_C

This differs from Eq 6, even if you assumed eps_T = 0.

Eq 15: Should this be epsilon_c rather than ∆c? Otherwise, define ∆c? 

Eq 15: The R << 1 case is clear.  It’s not clear to me how you are getting the R >> 1 limit from 13.  It seems like it should just be - (tau^-1 eps_c ) / R, since for limiting cases Rj is approximately equal to R0 (?), and thus the denominator reduces to 1 + (j^4 + 1) R + R + R^2, which would be dominated by R^2 for large R. Even including the j^4 R in the denominator doesn’t get me to the expression provided in Eq 15. 

284: You mean consistent with a specific real inversion? There are many inversion systems beyond Varon 2023a, for different observational datasets, trace-gases, resolutions, and transport models. It’s not clear to me that R = 0.2 is common across all such inversions. 

387: Outperforms in terms of error, but also probably is more expensive, in terms of computing K?  So there’s a tradeoff here.

382: It seems like p should change with the number of buffer grid cells.  With only 10 buffer clusters, it should be much larger, compared to the gridded buffer? 

Minor edits: 

16: estimates —> estimate

22: “grid cell fluxes outside the main of interest” isn’t very clear — for a regional inversion, wouldn’t the domain of interest be the entire region, and anything outside the BCs?  I suppose this distinction will become clear upon reading the paper, but at present for a reader skimming the abstract alone, it could be clearer.  The abstract is far from lengthy. 

27: “background concentrations given by boundary conditions” is a bit of an odd expression. Perhaps more direct to say “assuming boundary conditions”  

30: What is a “boundary condition concentration”?  I guess this is just a bit of imprecision in the language used here. I’d assume a boundary condition is a lateral flux into the domain.  This depends upon some assumed concentrations at the boundary, which would be “boundary concentrations”, not boundary condition concentrations? BC biases could be owing to the transport as well as the boundary concentration, yes? Some care with the wording regarding this throughout could be useful.  For example, line 34 would b e “Boundary concentrations used to define BCs can be provided by…”

40: And the magnitude of the inflow itself, so perhaps better to say “inflow conditions.”

34-36: Another common technique is to define background concentrations via statistical analysis of the lowest concentrations observed by the measurements used for the inversion.  

91: (DOFS), i.e., the number of…

98: as the diagnostic that estimates

102: It’s still not particularly clear what you mean by these two different methods, but maybe that will come shortly. 

123: the change in the influence —> the influence 

173: What equations are used to calculate the error for the intermediate case? 

196:  p..

198: Where does Eq 10 come from? 

391: Is this using Eq 9 , 13 or 15? 

S36: You mean Eq 9?