Publisher’s note: this comment was edited on 6 January 2023. The following text is not identical to the original comment, but the adjustments were minor without effect on the scientific meaning.

Review of “High mid-Holocene accumulation rates over West Antarctica inferred from a pervasive ice-penetrating radar reflector” by Bodart et al

Summary and main concerns:

The authors present an analysis of a ubiquitous internal reflecting horizon across much of West Antarctica dated to 4.72 ka, which they use along with a simple 1-D ice flow model to infer a spatially varying, time-averaged accumulation rate for 4.72 ka to present. They compare this estimate with model- and ice core-based estimates for present-day accumulation rates and find that the 4.72ka–present accumulation rate was roughly 20% higher than more recent accumulation rates, with some interesting spatial details. Most notably, while the overall 4.72–present accumulation rates were higher than present day, the sites closest to the coast (at Thwaites and Pine Island glaciers) show the opposite pattern.

Overall I think this is a good manuscript and will be suitable for publication in The Cryosphere. I have listed a number of comments on details of the figures and text that should be addressed in the section below. The biggest weakness of the paper is the lack of an attempt to quantify uncertainty due to the very simple flow model that is used here. While assimilating these data into a 2- or 3-D model with higher-order dynamics is obviously far beyond the scope of this study, the Local Layer Approximation and the 1-D Nye model come with many assumptions and simplifications that I think are likely to yield larger uncertainties than the authors assume. However, the benefit of such a simple model is that it is presumably inexpensive and relatively fast to run. Thus, a metric of model uncertainty could potentially be calculated by sub-sampling the dataset and repeating the analysis many times. If formally quantifying model uncertainty is infeasible, it is important to at least demonstrate that the results are insensitive to the choice of upper bounds on D, normalized depth, and horizontal strain rate.

Waddington et al. (2007) define the LLA as being suitable when D << 1, rather than just D < 1 as is used here. Because of the imprecision of this definition, it is not clear that even D = 0.25 is sufficiently small for the LLA to hold. MacGregor et al. (2009) seem to think even D = 0.1 is too high to satisfy the D << 1 criterion: “Although we are using the three shallowest spatially extensive internal layers in the radar data, nearly all of our study area (96–98%) has D values >0.1 that do not satisfy the D << 1 criterion, suggesting that the LLA may not accurately infer accumulation rates.” Their results also suggest that D might not actually be a good metric for determining where the LLA is suitable: “The values of D shown in Figures 3c and 4b suggest that the LLA is generally not suitable for flowband 1, but the small relative difference between b_LLA and b_fb for flowband 1 (5%) suggests that the LLA is acceptable.” They also note that while the LLA does an acceptable job in some circumstances, a flow-band-based  inverse model yields better results that differ significantly from the LLA in most cases. Thus, the validation of the LLA approach in this paper needs to be extended significantly beyond what is currently presented.

I am also skeptical that the LLA will be accurate in the lower regions of Thwaites and PIG within the model domain, where velocities exceed 300 m/yr, or in the regions where the 4.72 ka layer is deep in the ice column. It is also notable that the largest absolute values of Δb seem to occur in or near these areas. Care should be taken to show that there is not a significant correlation between Δb and D (or velocities or horizontal strain rates), and to evaluate whether the LLA applies when the layer is deep.

Specific comments:

L 71: Important to note that only some modeling studies do this for the time-scales relevant here. For instance, studies using higher-order physics and high resolution cannot hope to use LGM-to-present reconstructions to calibrate and validate their models and instead rely on usually just a few years to decades of observations. And it’s certainly fair to question whether lower-order ice sheet models like those used by the studies cited here should even be used to make predictions of WAIS changes over shorter (centennial) timescales.

Figure 1 needs an inset context map showing location of this area within Antarctica

Section 2.2.1 and 2.2.2: Out of curiosity, could this model be used to place bounds on accumulation rates where D≥1?

L 198: This is not correct. Waddington et al (2007) state that the LLA is valid only where D << 1, not where D ≤ 1 (see text below their eq 32). MacGregor et al. (2009) follow this definition. In contrast, MacGregor et al. (2016) do use D=1 as their threshold for the LLA, although I don’t see where they explain this choice.

Use of > and < on Figure S2 color bar limits are inconsistent

Use of sigma for strain rates is unusual. Can you use epsilon_zz for vertical and epsilon_xx for longitudinal strain rates?

Figure 2b shows that at some places the layer is ~80% of the ice depth. This goes against the assumptions in L 215–221: “where we can be reasonably confident that the ice sheet

has remained close to steady-state and where IRHs are likely shallow enough not to have sustained appreciable disturbances that would affect the Nye model assumptions “. It doesn’t look like the analysis is limited to depths shallower than 40% of the ice thickness. See other comments above and below about how the LLA needs more rigorous validation.

L223–226: But the analysis is extended to these areas of faster flow anyway, right? I don’t see any evidence that these areas are left out of the analysis.

Figure 3: Wording is a bit ambiguous in L306–307. Make it explicit that ice core accumulation rates are time-averaged for each core, not averaged across all 79 cores.

Showing the difference w/ cores in Figures 3 and 4 is helpful, but a lot of the information is obscured where the cores are very close together. It might be helpful to visualize this as a scatter plot as well as a map view. It also looks like there are a fair number of cores that give the opposite sign to the data in the map, but it’s not possible to tell from these figures how prevalent that is.

Figure 3c: inconsistent use of > and < between top and bottom of this color bar

L 342–345: This could also be an effect of the LLA being less appropriate at these lower elevations, where ice flow just happens to also be faster.

Really interesting that 4.7ka–present accumulation rates on Thwaites and PIG are so much lower than modern! If that’s real, presumably that’s an elevation effect. It doesn’t seem like this is discussed in much detail, and you focus more on the WAIS Divide story. However, this is potentially a more interesting result than the result for WAIS Divide, so I’d encourage you to add more about this to the Discussion.

It seems like the biggest weakness of this study is the inability to quantify uncertainty due to the Nye model. Is there a way to estimate that? For instance, there might be areas within the domain where the Nye model actually doesn’t work well, like in the tail of the normalized depth distribution in Fig 2c, or for other unknown reasons. Is there a way to re-do this analysis many times using sub-samples of the dataset to quantify the uncertainty due to the model? This could possibly be achieved by random sub-sampling, or by testing different thresholds of D or normalized Z. It would also be instructive to plot the inferred accumulation rate against D to ensure that D does not explain a significant amount of the variance in accumulation.

Can you report p-values for the comparison w/ RACMO as well as the comparison w/ ice cores?

L 323–325: I don’t understand the distinction between the two IRH-inferred accumulation rates in this sentence. Is the first one using a bilinear interpolation or some area average over the WD14 site?

Figure 5: As you allude to around L 347, I think it would be helpful to break this down into the IMIS and THW-PIG constituents. I think you should leave the current curves on the figure, but also add curves representing the same analysis for just IMIS and THW-PIG, respectively.

L 353: What does the “Fb” superscript here indicate?

L 373: This gives confidence that the Nye model works well at the ice divide, but that’s where it is most likely to work given its assumptions. This does not necessarily mean that it works well across the whole domain, where flow starts to deviate from divide flow.

L 376: “This also suggests that the WD14 Ice Core suitably represents atmospheric conditions across the wider WD.” How so? Can you explain further, and refer to a figure or table to help the reader understand the rationale?

L 404: Your comparison is between mid- and late-Holocene/modern. There is nothing here to suggest that the accumulation rate at 4.72 ka was higher than >4.72 ka.

L 420–422: Once again, I don’t see how the work shown here validates the WD14 record as a proxy for West Antarctica, and this is the first place in which the temporal variability of the WD14 record is brought up. This needs to be re-evaluated or explained more clearly.

L 445–450: Similar to comment ~L404: this argument is pointing the wrong way in time. Your inferred accumulation rates are higher than modern, but that does not imply that they represented an increase in accumulation rates relative to the earlier Holocene, which is what would be relevant for ice dynamics.

L 455: There is a recent TCD preprint that does suggest moderate grounding line readvance in the late Holocene: https://tc.copernicus.org/preprints/tc-2022-172/

L 465: It would be almost impossible for there to be lower accumulation rates at the coast than at the interior in Antarctica in the absence of >0°C air temperatures.

L 471–472: Same issue as above. How do your results show this?

The argument that ice sheet models need to account for time-varying accumulation rates seems like a straw man. Most ice sheet models do this in one way or another. While there are plenty of details about implementation that could be discussed, the point is not really brought up strongly elsewhere and feels out of place being stated so strongly in the conclusion. It has long been known that accumulation rates have changed from the LGM to present, and while the work presented here is a very useful and thorough addition to that story, it is not qualitatively different from what was known. While I agree that time-varying SMB should be taken into account, most models examining this period already do this. The difference between their various methods and the 18% difference between modern and average late-Holocene SMB is probably less important than model representation of ocean temperatures, sub-shelf melt rates, glacial isostatic adjustment, basal friction, low resolution and lower-order physics necessary for multi-millennium simulations, and other poorly understood processes that add enormous uncertainties to ice sheet modeling. It’s certainly worth discussing, but I would recommend de-emphasizing this aspect as a primary conclusion of the paper, since no analysis regarding the impact of the change in accumulation rates on model results is presented here. What would be more helpful is to emphasize that since the WD14 record does seem to be representative of the region of interest, SMB forcing for models can potentially use a modern SMB spatial pattern (e.g., from RACMO) that is modulated over time according to the WD14 data.

Supplement:

paragraph above eq S1: Is “decimated” the correct word? What interpolation techniques were used to sample from native resolution to 1km, and from 1km back to 5km?

Fig S1 caption: Last sentence should say “where D<1”, correct?