This is a clearly written paper with a controversial conclusion. The fact that temporal trends in the two polar regions are significantly different challenge our undrestanding on eithe ratmospheric process or the presented firn modeling efforts. This is adequately discussed in the manuscript. Hopefully, the paper will stimulate further reserach to solve this unexpected result. Presentation quality, language and level of scientific argumentation are excellent. The paper should be published in ACP

Review of Patterson et al. 2023, Reconstructing atmospheric H2 over the past century from bi-polar firn air records

Patterson and co-authors present new firn air-based atmospheric reconstructions of atmospheric H2 over the last century, combining data from several sites and both hemispheres. The difficulty in reconstructing atmospheric H2 from firn and ice core records is the molecular size-dependent bubble close-off fractionation. This fractionation results in H2 depleting in mature ice, and a corresponding enrichment in the firn air. The authors present several ways of correcting for this firn artifact.

I find their reconstructions convincing and worthy of publication in Climate of the Past. The work seems technically sound, and has been described clearly. I have a few suggestions I would like the authors to consider, in particular regarding the structure of the manuscript.

(1), I wonder why the authors chose to present the “regular” UCI_2 inversions as their main result (Figs 1-4), and then the various scenarios with altered pore close-off schemes as some sort of alternate or special case (Figs. 5-9). In my mind, one cannot reconstruct atmospheric H2 meaningfully without getting the close-off size fractionation correctly, and therefore I think the “altered” scenarios are far more convincing. Case in point, these “altered” scenarios provide a good fit to the atmospheric flask observations that are the gold standard in atmospheric trace gas reconstruction.

If this were my paper, I would not bother with the regular UCI_2 model, as it is clear it cannot fit the Neon data nor the direct atmospheric flask data. It is customary to calibrate (or tune) the diffusivity profile to trace gas records; in this situation I would consider fitting the Neon data as part of the model calibration (tuning). I would have only presented the “altered” scenarios. This would also result in a much shorter and more focused manuscript. I suspect the authors trust the altered scenarios better themselves, as this is what they plot in their last summary figure (Fig. 9). Perhaps this review can be a justification for the authors to make this simplification of the manuscript. On the other hand, I realize this is a big revision I am suggesting, and I would still be supportive of publishing this paper if the authors decide to keep the current structure.

At the very least the authors should communicate more clearly which inversion they believe to be most realisitic, so that users of these reconstructions know which one to plot. I would advocate strongly that the most realistic ones are those that can fit the dNe/N2.

(2) The authors apply a gravitational fractionation correction to the data via the d15N data. It would be trivial to similarly apply a close-off fractionation to the data via the dNe/N2 data (assuming that the dH2/N2 fractionation is the same as the dNe/N2, which the authors assume already). In this way the bubble close-off would be included empirically, and would not have to be modeled. Can you add an inversion using such an empirical data correction? My prediction would be that it matches the “altered” solutions.

(3) The discussion around the close off fractionation (section 7) is not as clear as it could have been. A few suggestions:

(3a) Can you add a panel to Fig. 7 showing the bubble pressure in the model? That is ultimately what drives the fractionation, so it is critical to have this information.

(3b) I don’t understand why the “reduced compression” scenario would result in more dNe/N2 fractionation. That makes no sense physically, as halving the pressure should also half the dNe/N2 anomaly. From Fig. 7a I suspect that the authors also altered the closed porosity parameterization, and that this is what drives the enhanced dNe/N2. Please check/confirm.

(3c) I am not surprised that the regular NEEM closed porosity parameterization gives a poor result for dNe/N2, given that it closes off much too deep (resulting in less pressurization). The Goujon/Martinerie close-off density at that site was artificially increased to match the field observation of the deepest pumping depth. It also ignores layering; including layering automatically results in some shallower trapping that will increase the dNe/N2 anomaly (as shown in Fig. 7a). Have you tried using the porosity parameterization from Mitchell et al. (2015), with the suggested close-off density from that paper? That may solve some of these problems. That parameterization does not produce an abrupt full bubble closure as may be required in some model architectures – this can be added manually perhaps.

(3d) The solutions from all three “altered” scenarios are virtually identical. Could this be because they start trapping bubbles at the same depth (Fig. 7)?

(3e) Could you extend the plot in Fig. 7 further down, to for example 90 m? Currently we cannot evaluate how sudden or deep the full bubble close-off occurs (the point where the closed porosity starts to decrease). Sometimes I find it more useful to plot the closed pore fraction, rather than the close porosity itself.

(4) Can you elaborate on using a Green’s function approach in the presence of bubble close off fractionation? Green’s functions assume a linear system response (the sum of two solutions is also a solution to the diffusion equation). Is this true in the presence of close-off fractionation? I suspect it is, but I am not entirely sure. Is the area of under the Green’s function greater than 1 in this case?

Other comments:

L33: “second-most”

L41: consider replacing “modern” with “present-day”

L45: Do you have a reference for the OH radical?

L52: Consider also adding a reference to Solomon et al. 2010, who first clearly describe greenhouse forcing from stratospheric H2O

L62: These trends are not very robust, and rely on single year anomalies. Also, how well are Khalil and Rasmussen calibrated with NOAA/GML?

L79-84: What about Antarctica? Those SH reconstructions are treated as somewhat of an afterthought in this paper, despite the topic being bipolar H2.

Section 2.2: Normally the diffusivity profile is somewhat model-dependent. Can you simply apply the profile from a different model?

L134: same depth “were” averaged

L136: intense seasonality: is it much deeper for H2 than for other gases like CO2?

L144: different sampling dates?

L154: This upper 5m is often called the convective zone

L154: this is often called wind pumping or just ventilation. Convection happens in winter when the surface is cold.

L158: that is a very small time step! I typically run my firn air model with a timestep of one week or so. Is this needed to keep the forward Euler scheme stable?

L170: Which parameterization? Goujon? Schwander? Mitchell?

L175: The d15N data also have thermal fractionation in them. How do you deal with this?

L189: Confirm that the model is coded in volumetric concentrations, rather than in ppm. Most models work in ppm, I believe.

L197: gas phase diffusivity: do you mean free air diffusivity?

L202: This is not really eddy mixing of course, though mathematically it is similar. This is more correctly described as dispersive mixing (Buizert and Severinghaus, 2016)

L201-206: so the advection is coded differently in the diffusive and lock-in zones? Is it a velocity term in the former, and a box-shuffling scheme in the latter? Does this conserve mass at the boundary? Do you account for the fact that there is backflow in the lock-in zone due to compaction?

Equations 4-6: How is this implemented? The x_n and P_bubble terms occur in all three equations, so you cannot simply solve them. Is this done iteratively? Or is there a typo in the equations?

L235: parenthesis ) missing after Rommelaere citation

Equation 7: Can you add the arguments to the variables here to clarify? For example, G(z,t) etc

Eq 8: What is N? Normal distribution? What is the time step i? 1 year?

L267: There is no artificial smoothing, but instead the parameter beta in the autoregression. Isn’t that just the same with a different name?

L275: With such cut-off depths, you have only 1 to 8 m in the diffusive zones. Can you confirm?

L 275: Can one instead add a seasonal cycle to the atmospheric history/ inversion?

Figures 1 and 2: is it possible to plot the firn data on the plot against their mean age?

L354: Can you add the Alert data to Fig. 3?

L357: Do you have a reference there for the ENSO connection? What about 1989?

L397: suggested that “the” maximum…

L409-410: But the maximum mostly disappears when you account for the close-off fractionation

Fig. 4: is it possible to plot the firn air data against their mean (or even effective) ages on the figure? Possibly empirically corrected for close-off fractionation? I always find this extremely helpful, as it allows the reader to visualize the data density and the degree of smoothing. As an example of plotting in this style, look at the recent Ghosh et al. (2023, https://doi.org/10.1029/2022JD038281)

L468: Have you also compared to Mitchell et al. 2015? The Goujon and Schwander parameterizations do not account for density layering, and are commonly applied incorrectly (i.e., they were derived on cm-scale hand-samples and are applied to m-scale bulk density; because of density layering this is technically incorrect. Layering broadens the depth of bubble trapping).

Fig. 5: It seems clear that the UCI_2 model does not build up pressure fast enough in the pores to expel fugitive gases. Shallower trapping seems like the obvious explanation to me, and tuning the closed-porosity parameterization makes sense to me as the first strategy. The Mitchell et al. 2015 parameterizations has two parameters that can be tuned.

Fig. 5: can you also show the fit to the Antarctic dNe/N2 data somewhere? These sites are used, and so the reader will wonder whether the modeling can fit those data.

L520+: can you give a plot of the s_c parameter? This is hard to visualize.

Eq 13: you have different numerical advection schemes in the diffusion and lock-in zones, right? How does this impact the implementation of the bubble pressurization?

L547: You seek to explain the effect via the bubble compression rate only, but then alter the porosity profile after all to get something similar to the first scenario. It is the porosity tuning that makes you fit the data, not the reduced bubble compression rate – if anything the latter should make it harder to fit the data. I think this is confusing to the reader. Why not instead conclude that solely altering the pressurization rate (which this scenario ostensibly represents) does not improve the fit?

Fig. 7: can you add  plots of the NEEM and Summit closed pore pressures? That seems needed to evaluate the fugitive gas enrichment.

L678: what about H2 artifacts during flask storage?