Comment on tc-2021-394 for the role of anthropogenic climate change in glacier

The overarching goal of this study is to provide a framework to attribute glacier retreat to anthropogenic climate change. The authors seek to do this by performing ensemble simulations of glacier retreat in idealized geometries driven by quasi-random climate variability. Overall, I think the study is very well written and illustrated. The figures were easy to read and interpret and the text provided sufficient motivation and narrative structure to follow the thread. I have some overarching comments, some boring if highly technical comments and some minor comments about wording, but overall my comments are minor.

past this implies glaciers should all be in their retreated position(s) and we wouldn't observe any glaciers in their more advanced positions. This isn't really a problem for the study because the climate, has not been stationary stochastic and has variability on a range of time scales. This raises two questions.
1. As the authors point out, modern glaciers are responding to both past and present climate forcing. We know that glaciers have advanced during colder periods. For example, some glaciers may have advanced during the Little Ice Age and then retreated and, depending on glacier size, glacier response to the Little Ice Age and other climate anomalies would overlap with the anthropogenic climate interval. However, the model clearly shows rapid retreat with little advance. This raises the question of whether the model (or perhaps geometry?) is capable of simulating prior glacier advance. If advance is not possible, then it seems possible that the model is overestimating the probability of glacier retreat in response to climate forcing (?), potentially biasing the statistical inference. To put this another way, in the authors model the glaciers will eventually retreat irrespective of anthropogenic climate forcing and the only thing that warming does is increase the probability that this occurs sooner. In this scenario, anthropogenic climate change only affects rates and not states (i.e., the time of retreat can be accelerated by warming, but retreat is ultimately going to happen irrespective). It would be more satisfying intellectually if the authors could turn around and also attribute glacier advance to periods when the climate was colder, like we have observed in the historical and paleorecords.
2. The authors break the probability distributions into a component related to (random) natural variability of the climate forcing and a component related to parameter uncertainty. This is fairly standard, at least in the glaciological literature and it follows from numerous studies in engineering. However, it makes a potentially large assumption: that we understand the system well enough that the model uncertainty largely derives from a handful of parameters that are imperfectly known. There is another possibility that also has to be considered which is that the underlying parameterizations are either not complete or fail in different climate scenarios. This seems especially relevant when dealing with submarine melt, iceberg calving, shear margin weakening, subglacial hydrology, etc, none of which are especially well understood. To be clear, my understanding is that the entire formalism presented here can be applied to any model irrespective of the models fidelity. For example, ca 2000 one could apply this same method using Shallow Ice Approximation models that don't account for longitudinal stresses or marine ice sheet instability. These models would require much more oomph from the climate variability to drive retreat because they lack crucial physics. But the same formalism would allow "attribution". Hence, a crucial point made by Shepherd (2021) is that we also have to consider all the alternative hypotheses that could also account for the observations. Shepherd (2021) described how to do this using Bayesian analysis through the use of the "complement". The trick is that one can formally include how much confidence we have in the model vs alternative models/explanations. I don't propose that the authors utilize this approach here, but I would like to make sure that they are aware of it and urge them to consider the possibility that their model might not be as physically robust as one might assume from the discussion in the text. I will note that this is gently hinted at near line 215, but it does seem important to emphasize that the attribution is very sensitive to model assumptions and ultimately, this effect needs to be quantified.

Technical comments:
1. How do the authors define noise and what does it mean for the forcing to bè`r andom''? My understanding is that the authors assume mass balance has a secular component with zero-mean fluctuations super-imposed. But I'm not entirely sure how the fluctuations are defined and I would encourage the authors to add additional details and equations about how the noise is created in the supplementary materials. More concretely, I take it that noise is added to the surface mass balance? For a zero-mean Gaussian process, the noise is not smooth and differentiable, so we would then need to integrate it in the form: $h(t+\Delta t) = h(t) + \Delta t (f(t)+S(t)) + \sqrt(\Delta t) \sigma(t,h)$ where \sigma(t,h) is the standard deviation of the Gaussian process and S(t) is the secular component and I have defined f(t) as the divergence term in the mass balance (or other terms in the equation). Note that the random noise term is multiplied by the square root of the time step in a Brownian process. There is a literature on integrating stochastic differential equations using colored (as opposed to white) noise, but that far exceeds my mathematical acumen. It would be helpful to me to see more details summarizing how the noise is created and how the stochastic differential equation is then integrated along with demonstrations of numerical convergence using both varying time step size and grid resolution. I don't request a host of convergence studies added to the paper, but a few sentences explaining that they were done and the results of the convergence experiences.
2. The frontal ablation parameterization is intriguing, but I have some questions and comments about this. As I understand it, this approach involves applying a large, negative surface mass balance localized at the last grid point at the grounding line and labeling this a "flux" or frontal ablation term. This seems intuitive at first: the flux term is removing ice at the terminus. The large frontal ablation causes a surface slope between the last two grid points in the (discretized) model. That this is in fact a surface ablation parameter can be seen by moving the frontal ablation to the right hand side of the ice thickness equation. For example, writing the flux as q, an upwind finite difference has the form: $h(x,t+\Delta t) = h(x,t) + \frac{q(x) -q(x-\Delta x)}{\Delta x}\Delta t + \frac{h \dot m}{\Delta x}\Delta t + S(t)\Delta t$ Note that the second to last term is the frontal ablation term. In the limit that $\Delta x$ becomes small, the surface ablation term becomes large, leading to an increased effective surface mass balance at that point. I have messed around with this type of parameterization a lot in the past (e.g., Bassis et al., 2017) and could not convince it to converge numerically under grid refinement. Instead, this type of parameterization created an unphysical singularity in the slope/thickness of the glacier that became larger and larger as the grid spacing became finer and finer. To cure the singularity, I had to regularize the frontal ablation term, recognizing that the surface ablation needs to be spread out over a characteristic length scale (I used 1 ice thickness). Doing this cured the lack of convergence and provided more physical surface slopes when using small grid spacing. But the results will depend modestly on the regularization scheme. Because of my experience, I would recommend considering a numerical convergence study to assess if the results are independent of grid resolution, time stepping, etc. This is not to say that the authors scheme is problematic, but it would be reassuring to provide some additional tests. To be honest, the entire attribution framework would still work even if the model does depend on the grid spacing. It would just emphasize my previous point that a real attribution requires some estimate of our confidence in model physics and numerics.
Bassis, J., Petersen, S. & Mac Cathles, L. Heinrich events triggered by ocean forcing and modulated by isostatic adjustment.Nature 542, 332-334 (2017) 3. I think my biggest recommendation is that the authors conduct some numerical convergence studies. In my opinion, numerical convergence studies are like brushing your teeth: unglamorous, but essential hygiene that needs to be regularly performed to avoid unpleasant surprises. This is often done and then forgotten about. Please tell readers what you have done even if you don't show it.

Minor comments:
Line 275 and elsewhere: While->Although. While is technically supposed to refer to time.
Introduction: Grounding line retreat in WAIS maybe related to the MISI (although also ocean forcing!), which is tied to retrograde bed slope. But there are other types of retreat. For example, the disintegration of ice shelves in the Antarctic Peninsula is not tied to bed slope (because the ice shelves are freely floating). Similarly, retreat of Petermann Ice Tongue is also not tied to bed slope. I think the discussion here is mainly focused on Greenland and grounded glaciers. This might be something worth emphasizing.
Bed topography in Figure 1 looks like it is piecewise continuous, but not differentiable. This can create numerical issues and problems with numerical convergence in models that assume the ice thickness is smooth and differentiable.
Line 135: Out of curiosity, why not use a one sided probability distribution (e.g, lognormal) that naturally avoids unphysical adding mass to the terminus? Figure 2 makes a key point: In a system that is close to a system with an instability, retreat always occurs and the only question is how long it takes for retreat to initiate. I don't know that there is strong evidence for this type of behaviors for glaciers. This might be partially because the climate is not stationary stochastic over this type of time scale.
Equation (4) appears to be missing an ice thickness on the right hand side? Please check.