Review of Sailer et al.: Ice core site considerations from modeling CO2 and O2/N2 ratio diffusion in interior East Antarctica

This paper presents simulations of how much CO2 and O2/N2 signals may be attenuated by gas diffusion within ice sheets preserving 1.5 million years of climate history. Also, using both an ice and heat flow model and a gas diffusion model, the authors predict regions across a broad area from the South Pole to Dome A where million-year-old ice is likely to preserve atmospheric signals with higher amplitudes. This study provides valuable information for selecting drilling sites for the NSF COLDEX project. Moreover, their estimates could be validated in the future through ongoing oldest ice core projects, which could in turn help constrain the diffusion coefficients of gas molecules in ice—parameters that are otherwise extremely difficult to measure.

Overall, I consider this paper suitable for publication in Climate of the Past after minor revisions. However, I recommend that all the points raised in this review be carefully addressed before the manuscript is accepted.

General comments:

1. Gas diffusion in ice is controlled by both the gas concentration gradient and temperature. This paper presents a set of sensitivity experiments by varying parameters such as accumulation rate, ice thickness, surface temperature, and geothermal heat flux (GHF), and it is a nice aspect of the work that the authors investigated how each of these parameters affects CO2 and O2 diffusion in the ice. However, to facilitate more intuitive understanding of the results, I suggest the authors include the vertical temperature profiles corresponding to the ice thicknesses used in each experiment. This would also help the reader better interpret the results of the experiments that assume the presence of stagnant ice near the bed.
2. Regarding the captions of figures and tables: in many cases (notably Fig. 4, Fig. 5, Fig. 6, Fig. 7, Fig. 8, Fig. 9, Table 3, and Table 4), the first sentence is descriptive. It would be more appropriate to begin each caption by stating clearly what the figure or table is showing. Moreover, several captions include explanations of the methodology or interpretation of the results, which should instead be included in the main text. I recommend that the authors review and revise all figure and table captions accordingly.
3. There are some citations to works that are in preparation, submitted, or in review. In general, such works should not be cited, as it is not possible to verify whether the citation or the associated discussion is appropriate. In particular, citing works in preparation is highly inappropriate. While I leave the final decision to the Editor, I believe that only works that are publicly accessible should be cited.
   
   Line 49 and more: Young et al., submitted
   
   Line 96: Singh et al., in prep
   
   Line 360: Fudge et al., in preparation
   
   Line 364: Parrenin et al., in review (there is no info in the reference list)

More specific comments

Lines 60 and 63: Line breaks are unnecessary here.

Line 69: The sentence “Feedbacks with the decreased…” is unclear. Please clarify the meaning and provide a proper citation.

Line 70: The reference to Stolper et al. (2016) is inappropriate here. It is related to reconstructing atmospheric O2 histories, not to dating. You should instead cite Kawamura et al. (2007, https://doi.org/10.1038/nature06015), who first used O2/N2 for dating. The correct reference for AICC2023 is Bouchet et al. (2023).

Line 72: Please consider adding citations to Kawamura et al. (2007) and Fujita et al. (2009, https://doi.org/10.1029/2008JF001143) here. They also discuss the mechanism for the relationship between O2/N2 and local summer insolation.

Line 73: The phrase “this age” is unclear. In addition, a citation is needed to support the statement that “This is the most reliable method.” For example, Oyabu et al. (2022, https://doi.org/10.1016/j.quascirev.2022.107754) provided the first evidence that O2/N2-based chronology is highly reliable, by demonstrating that O2/N2-derived ages show no phase lag relative to insolation cycles within the estimated uncertainties, based on comparison with U–Th-dated speleothem δ18Ο.

Lines 90 - 95 and Figure 2: The colors in Figure 2 make it difficult to distinguish between 2000, 2500, and 3000 m cases. Also, since the spatial extent is unclear, I suggest including a map of the entire Antarctic continent to show the location of the zoomed-in region. Please also plot the location of Dome A.

Line 97: What accumulation rate value is used for the past 4.7 ka?

Line 114 and below – Model description: It is helpful to present the numerical values of the model parameters (e.g., constants used in the equations) in a table. This would improve clarity and allow readers to reference them more easily.

Line 121: The value “2 cm/yr” is this in water equivalent or ice equivalent? 

Line 147: The Bereiter et al. (2014) model originally comes from Ikeda-Fukazawa et al. (2005), from which all the key parameters are derived. Please cite this work here.

Line 150: This part also needs a citation to Ikeda-Fukazawa et al. (2005).

Line 154: The correct references are Bazin et al. (2013) and Veres et al. (2013), not 2014.

Line 154 – 156: In this section as well, the authors should cite the original paper that provides the permeabilities, rather than Bereiter et al. (2014). The manuscript states that only O2 was included in the simulation. However, it is unclear how δO2/N2 values were calculated—was N2 assumed to remain constant at its initial value? The difference in diffusivity between O2 and N2 is less than one order of magnitude, about a factor of three based on the values from Salamatin et al. (2001). I am not convinced that it is justified to neglect N2 in the simulations. The authors should provide a clear justification or, preferably, a quantitative assessment of the impact of including N2 on the δO2/N2 signal damping. The manuscript states that “including N2 would slightly enhance the signal damping presented in later sections,” but I could not find any subsequent section in which the effect of including N2 was actually evaluated or quantified.

Start from line 187: This paragraph needs further clarification. Was the ice and heat flow model run over 1.5 Myr to derive average annual layer thickness and temperature every 50 kyr? I understand that the CO2 signal was prescribed with 5/4 cycles and the O2 signal with 5/2 cycles in the first 50 kyr interval, and that the diffusion was simulated under conditions where the annual layer thickness gradually decreases and the temperature increases with depth. Why was a 50 kyr interval chosen?

Line 189: The sentence mentions “Figure 3c and d,” but it’s unclear why they are cited here.

Line219: I do not understand why Bazin et al. (2013) and Veres et al. (2013) are cited here. Does this mean that the values in Table 1 were taken from these references?

Line 221: Similarly, I do not understand why Uemura et al. (2018) and Buizert et al. (2021) are cited here.

Table 1: While parameters for EDC are described in Bereiter et al. (2014), please explain more clearly in the text how the values for Dome Fuji (accumulation rate, ice thickness, surface temperature, GHF, p) were determined, and cite appropriate sources.

Figure 4: The lines with SDR values below 0.1 (corresponding to 0 to ~0.8 Ma) are nearly indistinguishable. Please consider adding a zoomed-in inset for this range to improve readability.

Line 279: The effect of ice thickness is difficult to interpret. While the O2 permeability from Salamatin is more sensitive to temperature than the CO2 permeability, the absence of temperature profiles in the ice sheet makes it difficult to fully understand the results. Please prove a more detailed explanation.

Figure 8: Please consider improving the color scale. In panel (a), the contrast in accumulation rate is barely visible. In panel (b), why is the color scale reversed compared to Figures 2 and 9?

Line 341: Please indicate in the figure that where Dome A is.

Line 396: “Our results show that CO2 SDR for 1.5 Ma ice does not exceed 13 % in the grid-north Foothills and averages 5 % (Figure 8k).” This sentence should include something like “with 50 mW/m² GHF” to clarify the condition of the result shown in Figure 8k.

Line 400: Did you also test a 20 kyr periodic CO2 signal? If so, please consider including a figure showing the results, either in Figure 3 or elsewhere.

Section 4.3: The experiment is unclear. In Case 2, does the temperature profile correspond to the upper part (0–2430 m) of the 2700 m case? Similarly, does Case 3 use the profile from 3000 m ice thickness truncated at 2700 m? If so, Case 2 includes both the thinning of the ice and lower basal temperatures, while Case 3 isolates the temperature effect. Please confirm and clarify.

Line 437: The accumulation rate is stated as 2 cm/yr here, but Table 3 lists it as 3 cm/yr.

Line 439: From “In Case 2, a thinner deforming...” I suggest inserting a line break to improve readability.

Line 440: Case 2 involves not just thinning of the ice but also changes in temperature? Please explain this more clearly.

Line 452: Oyabu et al. (2021) demonstrated that the permeability coefficients proposed by Salamatin et al. (2000) reproduced the smoothing of the O2/N2 signal in the Dome Fuji core well. However, their simulations were conducted over a limited temperature range. The temperature dependence of Salamatin’s permeability is quite strong, with estimates indicating that the permeability increases by approximately one order of magnitude for every 10°C increase. In contrast, the “Fast set” proposed by Ikeda-Fukazawa et al. (2005) exhibits behavior that approaches the “Slow set” at higher temperatures, such as near the base of the ice sheet.

While the reliability of these permeability estimates remains uncertain, if one takes an optimistic view, it could be argued that the use of Salamatin’s coefficients results in a more conservative estimate of signal preservation near the base of the ice sheet. Future measurements on actual ice cores may help clarify the temperature dependence of these permeability coefficients. It may be worth mentioning these points in the discussion.

Line 518: Adam Auton (2024). Red Blue Colormap (https://www.mathworks.com/matlabcentral/fileexchange/25536-red-bluecolormap), MATLAB Central File Exchange. Retrieved November 18, 2023.

This is not cited in the text.

Figure 1: The manuscript mentions “AICC2023,” but it is unclear whether this refers to the use of the AICC2023 age scale for the ice core chronology.

Review of Sailer et al., Ice Core Site considerations from modelling CO₂ and O₂/N₂ ratio diffusion in interior East Antarctica.

Sailer et al., have performed a critical analysis of the potential preservation of gas signals at locations in Antarctica containing the world's oldest ice with implications for existing and future ice coring efforts. I really enjoyed reading this study and particularly appreciated that the authors provided their code (which I have tested and look forward to using in the future!). By and large, the authors employ a previously established methods (Bereiter et al., 2014) and existing data on gas diffusion (Ahn et al, 2008; Ikeda-Fukazawa et al., 2005, etc) but expand the scope of previous work to include new, unexplored regions of Antarctica. In this effort, they also expand the range of parameters tested by Berieter et al., 2014 to investigate the sensitivity of their model to different temperatures, accumulation rates, geothermal heat flux, ice column thickness, etc. – effectively all the key parameters one would want to see varied when looking for old ice sites.  I would support publication in the Climate of the Past after some revisions, mostly minor.   

First, I will raise a few bigger picture questions and then go line-by-line.

Can this analysis be reconciled with the observations that there is O₂/N₂ variability in 1.5-million-year-old ice at Allan Hills? The major take-home from both Bereiter et al., 2014 and this study is that we should have lost all the 20-kyr and shorter variability.  Yet Yan et al., 2021 report variations in O₂/N₂ from very old ice at Allan Hills. I’ll admit that you’ll be comparing apples-to-oranges when jumping from a well-resolved deep ice core to a jumbled-up blue ice site where the ice (at least in the present) is colder than the ice buried ~3,000 meters below East Antarctica.  None-the-less, I can’t reconcile in my mind how O₂/N₂ variations would remain if the diffusion rates are indeed so fast. In fact, I even tried a quick test with your model by making an isothermal ice column at -30C.  The SDR of O₂/N₂ @1.5 million years age (with about 10,000 years per meter) was coming out quite high at 0.86.

On a somewhat related note, from Sailer et al. it’s not as clear as in Berieter et al., 2014 if we will have lost longer periods variation. Once limitation of this study is they only consider CO₂ changes on the obliquity timescale and O₂/N₂ on the precessional timescale.  I understand why the authors have chosen to impose the two different timescales to illustrate potential gas preservation.  On the other hand, it makes for a somewhat convoluted apples-to-oranges comparison for the reader during some later stages in the paper.  I hesitate to call for major revision, so I will only loosely suggest to the authors that they consider running experiments with the same timescale of forcings for both gases. One could then derive a parameter that describe the additional smoothing of O2/N2 compared to CO2. From the gas world perspective, yes, we want to know the absolute smoothing.  But we’d also want to know if the ratio of smoothing between gases changes with under different conditions.  

How does this study differ from Berieter et al., 2014 (if it all) in approach? From my reading, I believe the temperature-depth-age models are virtually indistinguishable. However, this would be nice to be spelled out exactly. Particularly if there are a few differences.

A few more tests to confirm performance of the model would be useful.  I recommend showing at least two examples of how the model performs for both a low-accumulation site with little melting (Dome Fuji?) and a high-accumulation site with high melting (WAIS Divide)?  In both cases, I believe the borehole temperature and melt-rates (possibly inferred not measured directly) data should be available.  It would be sufficient to this only for the review documents or in a supplemental figure.  I don’t think it would be necessary for the main body.  That said, the introductory material is somewhat lacking in a real-world grounding that could be better illustrated for non-ice core specialist.  For example, around lines 55-60, it’s taken for granted that reader has a good grasp how thinning and temperature evolve with depths.  This could benefit from an illustrative figure that shows some real-world examples. For example, the way Berieter et al., 2014 introduces the problem with real and modelled examples in Figure 1 is quite useful.

Also, on the subject of model testing, I took my best crack at an apples-to-apples comparison with a similar model we use in-house at BAS and found very good agreement in the modelled temperature and age with depth. Suffice it to say, the version of the model I used is mostly educational purposes such as teaching “where to find old ice?” exercises so a full model inter-comparison isn’t needed and well beyond the scope of this review.  However, I noticed one discrepancy could call for a little bit more description of your model.  In the scenario with -60C surface temperature, *2 cm a-1 a* accumulation and  3,000-meter thick ice column, I didn’t see any basal melting until the geothermal heat flux tipped over 55 mW a-1.

*Important aside, is that water or ice equivalent accumulation?

Upon further investigating and the running the code for your model I saw that the temperature of the bedrock (possibly down to a few kilometers is modelled). This is more sophisticated than the model I used and probably the main reason for the difference. It would be nice to know a little bit more about this portion of the model, particularly as the areas you eventually rule out for old ice exploration appear sensitive to the presence or absence of melting.  

Also, I struggled a bit to understand the iterative approach to solve the basal melt rate.  I believe a bit more detail is warranted along with a brief review of the implications of this approach. For examples, could you calculate how much heat is “lost” via conduction up into the ice and how much is “lost” via latent heat?  Also, I assume the implication is that the latent heat is indeed “lost” from the system?  As in, it is implied that a thin layer of water is flowing away from the site?

Finally, to circle back to my original suggestion of doing some model-data comparison, the “proof would be in the pudding”.  So if the model does well at simulating the temperature and melt-rates at Dome Fuji and WAIS Divide (or any other sites of your choosing)  then the reader will be more confident in the approach.

All the maps could be made more accessible.  I suggest you provide some more geographic context as I was a little bit lost as to the extent of the “COLDEX Search Region”. I suggest adding an Antarctic-wide map inset to Figure with the region currently covered in Figure 1 highlighted.  Dome A and Dome C play important roles in the paper.  Please show their location and also an arrow pointing in their direction (also Vostok?).  The subplots in Figure 8 (as long as they are same extent as Figure 1) are okay but seem to be missing lines of latitude and longitude.  I would also appreciate some more points of reference here.

Line-by-line comments:

Abstract:  “foothills” is not yet precisely defined and comes across as quite a colloquial term to use in the abstract. I recommend using a more precise term, or using the spatial information you provide about location between Dome A and South Poler, or add on a line like “…roughly equidistant between Dome A and South Pole, a region we call the “Foot Hills” of the Gamburstev Mountains”

Line 35:  “However, this method is limited, providing different results based on the species and location and requiring several assumptions. Köhler (2023) suggests that some of these assumptions may be incorrect by comparing reconstructions to a carbon cycle model.”  It would be helpful for a reader unfamiliar with Kohler et al to mention some of these key assumptions.  Otherwise, it sounds like a quite a broad, unsupported swipe at the boron method.  

Lines 55-65.  There are some very short, two-sentence, paragraphs here. I think it could be restructured into one, possible with some bullet points.

Line 74: “This is the most reliable method for dating ice cores of this age.”  A fairer statement would be that O2/N2 is one of the key pillars of dating as in practice all available information is used (e.g. Bouchet et al., 2023)

Line 150. “Bubbly ice resides in the upper ~1000 m where diffusive smoothing is unimportant due to thicker annual layers and colder temperatures.”  I would rephrase.  As shown by experimental evidence that underlies these results (Ahn et al., 2008) diffusion dues occur in bubble ice. In fact, one could argue that those rates derived in Ahn and most subsequent work don’t apply for clathrate ice in question here. I would say something like “although diffusion occurs within bubbly ice, the time spent within the bubble phase is relatively short (e.g. 25,000 years at a site like EDC) compared to the timescales of interest here (e.g 1,500,000).”

Please provide some more introduction to “fast” and “slow” datasets.  Why the large discrepancy? I wouldn’t expect the authors to solve the problem but a figure like presented in Bereiter et al, 2014 would be helpful. I found myself going back and forth between Bereiter et al., and this study quite a lot.

Table 1:  I know it is mentioned somewhere the text, but it should be reiterated that you appear to be using some sort of glacial-interglacial average for the sites. Otherwise, -60C at EDC jumps out as the reader as strangely cold. 

Figure 4.  “Current ice core measurements cannot be used to estimate diffusion in older ice.”   I would disagree with such an unequivocal statement. One take home (I had) from Bereiter et al., 2014 is that millennial-scale and faster variability can be used to estimate diffusion. A narrower statement, like “current orbital-scale variations in gases cannot be used to estimate diffusion” would be more apt.

Line 405: Does TAC diffuse? This is an interesting point.  Actually, I believe it could be easily tested with your model as Uchida et al., 2011 provide an estimate of whole air diffusion.  They are shown in a figure in Bereiter et al., 2014.  

Overall, very nice work.  I’m looking forward to seeing the revision and then the paper published.  Also, I’m excited to use this model!

All the best,

Thomas Bauska

Uchida T, Miyamoto A, Shin’yama A, Hondoh T. Crystal growth of air hydrates over 720 ka in Dome Fuji (Antarctica) ice cores: microscopic observations of morphological changes below 2000 m depth. Journal of Glaciology. 2011;57(206):1017-1026. doi:10.3189/002214311798843296

Yuzhen Yan et al. Ice core evidence for atmospheric oxygen decline since the Mid-Pleistocene transition.Sci. Adv.7,eabj9341(2021).DOI:10.1126/sciadv.abj9341