Please see attached report.

This paper details cosmogenic ¹⁰Be and ¹⁴C measurements in bedrock sampled right at the margin of four glaciers across the western US and Canada. The authors apply a Monte Carlo forward modeling technique to determine the most likely timing of glacial advance beyond the modern extent during the Holocene. At face value, the results yield different Holocene glacial histories among the four glaciers, which the authors suggest is unexpected given the uniformity in Holocene climate forcing across the region. To explain this unexpected heterogeneity, the authors hypothesize that the bedrock sampled at each glacier contains information about a unique part of the Holocene ice-margin history because these glaciers are now responding unequally to modern climate change. They attribute the difference in the magnitude of response among these glaciers in part to differences in glacier response time.

I enjoyed learning about the methods employed to interpret the cosmogenic ¹⁰Be and ¹⁴C data, found the figures compelling, and think the proposed hypothesis is interesting. I do have several major questions/comments/concerns that I believe need to be addressed before this work can be published, which I’ve detailed below. Many of my comments can likely be resolved by tightening up the organization and language. I also think the manuscript would benefit from including more details about some of the methodology and including a more comprehensive literature review in some places. I think these issues can certainly be addressed and hope that the feedback I’ve provided here is helpful.

Thank you for the opportunity to review this paper and I look forward to seeing where it goes!

General Comments

Methods: I like the approach used here taking forward modeling with Monte Carlo to estimate the most probable exposure histories represented by the cosmogenic-nuclide concentrations at each of the four glaciers. However, I was left with several questions about the methods used at several points in the paper, in particular in the description of the forward modeling. A lot of the information I was looking for is actually in the paper, but I felt it was buried in table/figure captions or was obscured by confusing (at least to me) language.

Response time calculations, near L103: How are mean glacier thickness maximum mass loss calculated? I think it is mentioned in the footnote of Table 1 that these values come from other studies. It would be helpful to cite their source in the text when you discuss calculating response time. Also, I think the wording above Equation 1 is backwards? Should it say mean glacier thickness divided by maximum mass loss from the terminus?

Apparent exposure age calculations: There is quite a bit of information buried in the footnotes of Table 2 that might be useful to bring into the text in Section 3.1. In particular, I think information about when you conducted fieldwork and the production rate used in exposure age calculations would be useful to have in the text. In the Table 2 caption, you say that you used the “global” production rate- do you mean the primary production dataset of Borchers et al. (2016), which is the default in v3 of the online exposure age calculator you use? If so, this production rate dataset does not include all production rate calibration data globally, so it might be better to refer to it as “the primary production dataset of Borchers et al. (2016).” Someone corrected me on this recently, I thought I’d pass it along!

Monte Carlo forward modeling and interpretation: I realize that the methods applied here are described in detail in the Supplement to Vickers et al. (2020), but I think this paper could include a bit more summary information so that the approach taken is clearer. I was left with several questions:

• How are the burial histories constructed? Is it totally random whether a given timestep has exposure or burial? L208-209 say “the scenarios include a range natural variability intended to simulate natural glacier length oscillations”, but I don’t really understand how that is implemented.
• Is each sample allowed to experience a different amount of erosion? It seems like it based on the last sentence in Section 3.2, but I couldn’t tell based on how the method was described. This is important if erosion is used to explain scatter in the datasets.
• How were production rates calculated with depth? Basically, what attenuation length was used for spallation and how were muon production rates calculated? There’s no mention of production pathways at all, but this is important, especially given the interpretation that some of the higher 14C-10Be ratios at JIF Glacier result from deep glacial erosion.
• The authors describe how the probability of exposure was calculated in the caption of Figure 6, but I think it warrants mentioning in the text as well.
• A statement on L278 says, “The probability of bedrock exposure at each site exhibits a quick transition from exposure to burial. For simplicity, we interpret burial to begin at the mid-point of this transition and report the time in which glaciers advanced past modern positions as approximate values,” but in the Figure 6 caption it says, “We interpret the change in probability of exposure from greater than 50% to less than 50% to indicate roughly when the glacier advanced past its modern size and bedrock samples were buried”. Maybe I am missing something, but I think these are different approaches for interpreting the timing of burial (i.e., none of these probabilities have mid-points at 50%). Based on the reported numbers and looking at Figure 6, I think the wording in the figure caption is more accurate?
• L220: “Samples with 14C-10Be ratios above the production ratio for surface exposure are excluded from the Monte Carlo forward model because they are theoretically impossible, though the possibility of high erosion rates exhuming nuclides produced in the subsurface is explored in the discussion” – It is discussed later in the paper that these ratios can be achieved with deep erosion, so in that context they’re not theoretically impossible. Are there model runs that achieve these ratios? If so, do these samples need to be excluded? Or are the erosion rates explored not high enough to achieve these ratios? This comment goes back to my question above about how production rates are calculated because if simplifications were made to the muon production calculation maybe there wouldn’t be elevated 14C-10Be ratios in your model results?
• Related to the above comment, the paragraph starting on L349 could maybe use a more technical explanation of why 10Be and 14C production change differently with depth?
• Is it permissible for some samples at a site to be buried while others are not? It seems like not, which is fine, but maybe you could include somewhere the distance of the samples to the 2021 terminus to further justify the assumption that all samples at a site should have the same exposure history.

Treatment of Subglacial Erosion: I appreciate the thorough discussion of subglacial erosion as it is crucial to interpreting the exposure age chronologies. First, I do have a few questions related to the Conness Glacier results are interpreted/presented.

L285-286: I wouldn’t consider erosion rates of 0.5 mm/yr to be “exceptionally high.” Rather, I think this falls pretty squarely into the range (and maybe even slightly lower than average) expected for mountain glaciers at mid-latitudes? See Cook et al. (2020) titled “The empirical basis for modelling glacial erosion rates”, which should be cited here and elsewhere when discussing expected/typical erosion rates (and references therein). In the same paragraph, the authors state, “all three glaciers expected to be more erosive than Conness”. Why is this? That statement seems to motivate the “cutoff” value of 0.5 mm/yr as the maximum reasonable erosion rate, so I think it deserves more of an explanation.

Based on the figure in the supplement, it’s clear that there are an abundance of exposure histories that replicate the nuclide concentrations at Conness when erosion rates are allowed to be higher than 0.5 mm/yr (alluded to in the main text as well), which is why the probability of exposure is <50-60% through the Holocene. I’m just wondering if, rather than burying these results in the supplement, could they be added to Figure 6 as an alternative explanation for the Conness nuclide concentrations, which might make this easier to discuss in the main text as well and explain why it’s not the preferred interpretation. It could even be a dotted line to make it clear it’s not the preferred interpretation and add that reminder in the figure caption. I don't want this to be confused with over-interpretation of very low nuclide concentrations, so maybe it's not the right solution, but something along these lines may provide some help for the reader in understanding the preferred explanation for these low concentrations. 

L412-414: “Subglacial quarrying is a plausible mechanism for deep erosion, but quarrying is spatially heterogeneous (Woodard et al., 2019); the measured concentrations at Conness Glacier are remarkably consistent and clustered tightly around the origin of the isochron plot.” I agree with this statement. Something that could help the authors rule out this quarrying scenario (as well as the talus scenario below) could be to put numbers on the range of erosional depths needed (not rates, but total depth of erosion during late Holocene burial) to achieve the measured nuclide concentrations at Conness. This might help with visualization of what size blocks would need to have been removed to yield such low nuclide concentrations (likely large enough that this is not plausible). Also, if the numbers end up being relatively uniform across samples, it’s probably unlikely that there were six uniformly large blocks quarried and you just happened to sample at those sites. I think it would be pretty easy to calculate from the model results by multiplying the preferred erosion rate by duration of Late Holocene burial and making another histogram like those already shown in the supplement?

One more comment about subglacial erosion: These results give estimates of subglacial erosion rates at several glaciers across Western North America, which is an important finding because subglacial erosion is something that is difficult to measure because it takes place beneath glaciers. If the authors wanted, they could devote a little space in the discussion to highlighting this finding and placing it in the context of the literature (this is sort of done when stating why the erosion rates at Conness shouldn’t be >0.5 mm/yr, but I feel like could first be in a more positive light for the other three glaciers, which may bolster the argument for lower erosion rates at Conness).

Hypothesis about non-uniform burial durations

The major conclusion of the paper is that Holocene retreat histories were similar across sites (as the authors expected) but the four glaciers have experienced different magnitudes of modern retreat so that the glaciers that have experienced more modern retreat are now revealing bedrock that spent more of the Holocene covered by ice. This is summarized beginning on L440: “We find the most parsimonious explanation of the non-uniform bedrock burial durations to be that glaciers advanced roughly in concert across North America, as predicted, but have retreated non-uniformly in the industrial era.” This is taken one step further to say that “modern glaciers are no longer behaving in concert like they did over the Holocene and are now uniquely out of sync”. I think this is an interesting hypothesis but there were a few key things missing for me to completely follow the argument. In some places, it feels like there is some extrapolation beyond the dataset that is not necessary.

First, this argument seems to rely on the set up that all four glaciers are expected to have behaved more or less in concert throughout the Holocene because the “first-order climate controls are similar” (~L60) across the region. The uniformity of Holocene climate change across the Western US and Canada is again stated briefly in the first paragraph of Section 5.2, but for the most part the reader is expected to take these statements at face value. Given that these glaciers span 22º latitude, occupy different climate zones/mountain ranges, and are at different elevations, a slightly more detailed lit review could be given to this topic to establish for the reader that these glaciers should have experienced similar Holocene climate and therefore are expected to have behaved in concert.  

The interpretation that modern retreat is tapping into different parts of the Holocene exposure history is really interesting. My issue is the definitive wording surrounding this interpretation. Because this explanation requires that each sampling location yields information about a unique part of the Holocene related to how much the modern glacier has retreated, I don’t think you can definitively say that these glaciers did (or didn’t) behave synchronously during the Holocene, which is stated several times. It’s possible I’m missing something here, but I think the manuscript would benefit from toning down some of the more definitive language in Section 5.3.

Organization and writing style: I recognize that my role as a reviewer is to comment on the scientific rigor more than the writing, which is of course up to the authors as approved by the editor. I do feel, however, that this paper would benefit from another pass at organization and wording. For example, I felt the Results section jumped around a bit between burial durations/erosion rates as derived from the two-nuclide diagram and the apparent ¹⁰Be exposure ages. Could results be restructured to first introduce apparent exposure ages, since those are the simplest to describe, then introduce where they plot on two-nuclide diagram? Similarly, the discussion could be organized so that the exposure history of each glacier is established before the non-uniformity of the records is discussed? When discussing the exposure history at each glacier, it could help to start with Konakee and Mammoth because they require the simplest explanations, then go into JIF and Conness which require more complexity. I recognize that there’s a lot to cover in this paper and that there’s not one right way to organize it, but these are things that as a reader might have worked better for me. I also think a lot of my comments above may have been clarified if some of the wording had been more precise, which will come up in some of my minor comments below.

Specific comments:

L19 and elsewhere: “rapidly retreating to their sizes today” – not sure what size means here – length, volume, thickness? I assume lateral extent? This phrasing comes up throughout the paper and could be clarified.

L61: “Glacier length is also a robust record of climate” – does this mean summer temperature?

L40: “how” is vague, does this mean the degree to which glaciers fluctuated? Response of glaciers to past warming/cooling?

L92-94: The authors state that glaciers across region began retreating from LIA moraines 1880 CE, but don’t say how that is known (like, what did the authors in the cited study do, what is the evidence, is it really that all glaciers in the western US began retreating by 1880?). Could there be a bit more explanation? In general throughout the paper I feel like there are statements like this (usually with citations, which is great!), but as a reader who is not an expert on when glaciers in the western US retreated from their LIA moraines, I just was looking for a bit more info about the literature. I'd just keep an eye out for these! 

L105: “we caveat this response” – I’m not sure caveat can be used as a verb like this?

L111: “lowest-sloping ice surface” is confusing – maybe, surface slope of only 5º, or, shallowest/gentlest slope?

L111: units for response time. Years?

L111-112: “most impacted by absolute area loss from climate change, but least impacted by percentage relative to its Holocene maximum at 70% of its LIA area” – I’m not sure I understand this. Absolute area has changed the most since the end of the LIA but still occupies 70% of its LIA extent? The phrasing, “LIA area loss”, used throughout Section 2, is confusing to me. Does this mean the glacier in 2021 was 70% of its LIA footprint?

L 112-113: “2003 area of JIF glacier mapped in Figure 1 is likely similar to area of the glacier at ~2 ka prior to the latest Holocene and LIA advances” – Is this stated in Clague? If so, could use more of a (brief) summary of what they did and how they came to this conclusion.

L143: Table 1 says that response time couldn’t calculated for Conness because of a lack of suitable data (which I assume to be mean thickness and maximum mass loss?). Paragraph here says “RGI mapping is too different from glacier size today”. I don’t understand what that means and why that precludes calculating response time. Also, could add something like “given the relationship between surface slope and response time, likely has response time similar to Kokanee Glacier”?

L182: Several samples, or several samples at each glacier?

L183: “we assume that all samples [within a catchment, at each glacier] experienced…”?

L263-264: “the Conness glacier … four glaciers preseted here” – these sentences basically say the same thing. Could condense.

L297: “results at all four glaciers best explained by glacier expansion from the early to late Holocene” – continuous expansion from the early to late Holocene? I don’t think these data show whether that glaciers were expanding or retracting through the entire Holocene, but do show that three of the glaciers were generally smaller than today in the early-to-mid Holocene and larger than today in the late Holocene.

General comment about presentation of apparent exposure ages: how is the cutoff of “three highest concentration samples” determined? Also, is this truly highest concentrations or oldest ages? Possible it’s the same for these sites because elevations/shielding are likely similar?

L303: Could use more context about the increasing glaciogenic sediment flux in Konrad and Clark, 1998 – how long is the record? Is there glaciogenic sediment the entire time? Could give more space to this comparison, as it seems like these data could potentially help support the preferred interpretation.

L367: What does minimum estimate mean here? That JIF likely advanced earlier than 2 ka?

L376: “retreat events did not extend beyond modern position of the glacier” – should this be that the glacier was not smaller than today during those retreat episodes?

L389: What does it mean for “the low 14 C10 Be ratio of MG-02 to be an artifact of the simplifications made for the isochron plot”?

L454: I thought the preferred interpretation was that Conness glacier samples were buried the whole Holocene. How did it recede past its minimum size of the Holocene? Wouldn’t that imply exposure at some point?

L462-463: “Kokanee Glacier has retreated 46% from its LIA extent, an intermediate amount between Conness (14%) and JIF (70%) Glaciers.” This makes it seem like Conness retreated the least and JIF retreated the most compared to its LIA extent? I think it’s the opposite?

L458: What does “nearly the least relative to its Holocene advance history” mean?

L459-461: This is an important sentence but I’m not sure I understood it – what is “the further back into its Holocene advance history it resides today” – revealing rock that spent less of the Holocene covered? 

Table 1: Does it make sense to provide mean glacier thickness and mass loss at terminus? I realize that these are from other studies and are included in a supplement al table, but that supplementary table could probably be eliminated by including these values in Table 1.  

Technical corrections:

L23 and elsewhere: minor, but I believe Icefield is all one word?

L90: moraines plural?

L191: I think spatially should just be spatial?

L370: inferring should be inferred?