review of “A probabilistic framework for quantifying the role of anthropogenic climate change in marine-terminating glacier retreats” by Christian et al.

This study considers the question of how to quantify the impact on glacier retreat of externally-forced climatic trends, in the context of internal climate variability. The study uses large ensembles of simulations of a simple glacier model to illustrate several points. (The members of each ensemble each having a different realisation of random climatic variability.) One conclusion is that the influence of a climatic trend on a single glacier geometry may be quantified by comparing ensembles of simulations with and without the trend. Another conclusion is that the influence of a climatic trend may be quantified by considering a population of glaciers with different geometries, all under a common forcing. These conclusions may also be stacked, considering a population of glacier geometries under a population of trends. The results illustrate a variety of philosophical and physical considerations that must be addressed in order to make statements about whether contemporary ice sheet retreat is linked to human behaviour.

My assessment is that this is a very nice study that deals elegantly with an important and very challenging problem. The paper is extremely lucid and thought-provoking throughout, and I feel it has taught me a lot.  One could perhaps argue that the conclusions are not particularly surprising, but I think that is a very good thing – in my opinion all the best science leads the reader to an understanding that feels intuitive and ties everything together.  In this paper, the authors have bought some of the key principles of climate-change attribution to glaciology, and I think the glaciology community could benefit hugely from taking these on board.  Further, I feel that the strong nonlinearity of the glacier dynamics studied here means that this attribution work has features that are rarely found elsewhere in the climate change literature.

I offer a few comments below but I have no reservations about this paper.  I think the authors should be allowed to consider my comments below and respond to them however they wish, including ignoring them entirely.  I’m happy to read the responses but don’t require this prior to publication.

Comments (in order of appearance in the text)

Introduction: The introduction frequently mentions Antarctica, but after a while in the following methods section it became clear that the simulations here only explicitly deal with marine-terminating glaciers.  The ice streams of main interest in Antarctica all have ice shelves.  I feel the paper could be clearer on the extent to which the results apply to the case in which an ice shelf buttresses the glacier.  Perhaps the authors can state at the very beginning of the paper that they don’t study that case, but that the ideas are broadly applicable, and maybe return to that in the discussion.

Section 2: I was a little confused about what the authors mean by frontal ablation.  To me, ‘ablation’ means melting.  However, maybe this is meant to be a combined melting+calving rate? But then, the model does not represent any floating ice and the ice is calved wherever it goes afloat, at a rate equal to the flow speed.  I assume the ablation rate is applied to the face after it has been calved at floatation.  I guess the ice retreats happen because an increase in ablation displaces the terminus back a little, then the ice speeds up, thinning the ice, and causing further retreat of the floatation line. So, I am a bit confused.  I assume I am just showing my inexperience here as this is a well-established model, but perhaps the authors can discuss this a little in the paper as I think the mechanism is relevant.

Section 3:  In Figure 3a, why do 68 of the ‘natural control’ experiments retreat during the decorrelation period that precedes the time period shown in the figure (i.e. only 932 are left out of 1000 started)? Does this imply the period shown, in which no retreats occur, is not representative?  I also note that the simulations in panel b are retreating right from the start of the period shown, so presumably these plotted retreats are just a continuation of the population of retreats that occured within that ensemble’s decorrelation period? Do these matters imply that the conclusions drawn in this section are dependent upon the decorrelation period chosen? 

Section 3.1: When varying model parameters, the authors find that the probability of retreat is a function of the distance of the steady-state terminus from a bed peak.  To be clear, are the authors saying that they think this relationship is causal – i.e. a larger displacement of the terminus is required to reach the peak and trigger retreat – or just a correlation – i.e. the underlying physics of the glacier have been changed in such a way as to enhance instability?

Section 3.2: (line ~260) I didn’t quite follow the physics here. I naively feel that an increase in discharge could enhance the advection of thicker ice towards the bedrock high, hence advancing the terminus. Is it always the case that an increase in discharge is the crucial destabilising factor?  I guess it is actually an increase in ice divergence that is needed to thin the ice and retreat the terminus? By the way I am aware of other literature that draws relevant conclusions concerning the frequency response of ice streams to climatic perturbations (e.g. 10.1098/rspa.2012.0180, 10.1002/2017GL075745) though that is only in the ice shelf-buttressed case.

Section 4:  (line ~302) I didn’t quite click with the language that a background trend makes the positive anomalies more persistent.  Which is more important:  the slow thinning of the glacier that I assume accompanies the (quadratic?) time-integrated ablation anomaly, or the fact that any given positive ablation anomaly is larger?

Section 4 General: The paper discusses the difficulty of defining a retreat metric, which I fully sympathise with.  In this section, the paper determines the effect of a climatic trend on the probability of a retreat of a given distance within a fixed time frame – e.g. before the end of a 150 year run.  Under this approach the probability is variable and the time frame is fixed.  I wondered if the authors had also considered the inverse approach – asking what is the ’time-to-emergence’ of a fixed probability of retreat.  E.g. if we choose to be interested in a 50% probability of retreat, the authors could determine how long any given trend would take to induce such a probability (compared to a no-trend scenario). This would have the advantage that the outcome is not a function of the arbitrary duration considered (replacing that with the arbitrarily chosen probability).  I recognise that the approach currently taken may be more appropriate to historical attribution, and the time to emergence idea is usually used for projections.  My guiding principle here is that under ANY nonzero climatic trend, eventually ALL glaciers will have retreated.  So, to me, a time-to-emergence metric reflects that situation.

Section 4 General:  Is there a significance test that needs to be applied here? If we see a difference in retreat probability between ensembles with and without a trend, or with two different trends, surely we need to determine that difference is statistically robust?  I cannot immediately think what would be the appropriate test, but I assume it would tell us what ensemble sizes are needed to establish a given retreat-probability difference between two ensembles at a stated confidence level.

Section 5 General: This section assumes that all glaciers in the population have identical climatic forcing, which seems a little restrictive to me.  For example in Greenland, all glaciers experience similar atmospheric conditions and far-field ocean forcing, but that is quite different to saying they have the same SMB and frontal ablation rates, which are determined by very local features such as ice topography, fjord geometry etc.  I believe the logic assumes that if neighbouring glaciers have different retreat history, that can only be caused by terminus bed geometry, which I don’t believe is always the case.

Section 5 Figure 7: I was initially surprised that the ensemble in panel a has only one member that advances.  I believe this is telling me that there is a statistically significant internally-generated trend in the climatic forcing (towards retreat).  This means that the ensemble is ‘primed’ such that when the external trend is added in panel b, lots of glaciers retreat. This is useful for illustrative purposes, but it is not mentioned and I think the authors should be open about this situation.  They could potentially add an internally generated trend line to panel a.  They could add a red dot to panel c illustrating that this chosen realisation sits above the mean fraction retreating for 30 m/y (I assume).  Probably the best thing would be to select a different realisation that has zero internally generated trend.

Section 5 General: As with section 4, what statistical testing would be required to demonstrate that a population of glaciers was retreating under climatic forcing, relative to a population fluctuating with no climatic trend, at a given confidence level? 

Section 5 General: As a closely related point to the one above, I found myself wondering what is wrong with just asking what fraction of glaciers in an area have retreated in the real world.  If enough Greenland glaciers are monitored, over a long enough time period, any net retreat implies a climatic trend in forcing must be important, does it not?  Then the question becomes how many glaciers and how long a time period need to be monitored to provide a given confidence level.  This re-states my ‘time of emergence’ point above.  I can’t quite link this concept to the work in section 5, but I bet the authors can. (Plus, I bet enough Greenland monitoring data are available to provide a pretty high confidence level.) 

General:  Even if the existence of important climatically driven changes in a glacier can be established, that does not imply that the climatic changes are anthropogenic.

Conclusions Line 521: natural fluctuations in climatic forcing

Summary.

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.

Overarching comments.

One of the results of the manuscript is that random climate fluctuations will eventually cause glaciers to retreat.  The time it takes to do so depends on the magnitude of the imposed noise.   Of course if glaciers had experienced a stationary stochastic climate in the 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 paleo-records.    

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 be ``random’’?  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, log-normal) 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 (1) in the supplement: I think the exponent is (1/n)-1 and not 1/(n-1).  Please check.

Equation (4) appears to be missing an ice thickness on the right hand side?  Please check.