We would like to thank the reviewer for taking the time to review our manuscript, and we appreciate the positive attitude towards it as well as the helpful comments to further improve it. We reply below by highlighting the reviewer's comments in italics above our answers.

My main comment is the discussion of bumps and no or slow retreat of the GL while the authorsuse SSA - plus no subglacial hydrology, and the Budd type friction law with a uniform friction coeffi-cient. I understand the choices for the experiment setup but would wish for a more detailed discussion how a higher order or full stokes simulation would affect the results on bumps and overdeepenings (e.g.SSA critised for a similar setup in Favier et al. (2012)). Furthermore, this concern is only raisedin sub-section 4.2. Before (Lines 217ff, 263, 367ff, . . . ), it is only mentioned that the reason for no retreat at bumps is the shallowness of the fjord.

We acknowledge that the SSA is not a full representation of the complex stress regime near thegrounding line. For weak beds and fast-flowing outlet glaciers, as we aim to mimic in our synthetic setup, the SSA is however a reasonable approximation and widely used in the glaciological literature. In the study cited above (Favier et al., 2012), we do not see any criticism for using the SSA in a similar setup. The comparison of their results using a full-Stokes model to a setup using the SSA refers to a study by Goldberg et al. (2009), but, as is stated in Favier et al. (2012): ”The authors [Goldberg et al. (2009)] also observed a merge between both grounded areas with advance rates of the grounding line situated around 70 m a−1, whereas in our case the rates are lower and between 15 and 30 m a−1. This discrepancy is not surprising since both experiments were performed under different configurations and are hence not directly comparable. For example, the SSA experiments were conducted with a null velocity prescribed along the lateral walls, whereas we imposed free slip.” (Favier et al., 2012). We do not see how this statement, or any other made in Favier et al. (2012), discount the appropriateness of using the SSA in a setup similar to ours. On the contrary, we find that the results presented there using a full-Stokes model agree favorably well with our SSA results for bumps and overdeepenings. In Fig. 11 of Favier et al. (2012), longitudinal deviatoric stresses are plotted alongside along-flow profiles through the glacier. The stress distributions show the same characteristics as in our results (Fig. 5), with compression on the stoss side of a bed rise, and extension on the lee side. Furthermore, the along-flow profile show that on top of the pinning point in Favier et al. (2012), the glacier surface shows a bulge which is exactly the same feature that we observe in our setup (Fig. 3c.). Therefore, we believe that the study by Favier et al. (2012) rather supports our methodology in that it shows that a full-Stokes or higher-order model would likely yield similar results as the setup with SSA used here, and that the conclusions would remain.

With this said, we will expand on the discussion on the influence of the ice-flow approximation, following the discussion above. This will be done in 4.2 Study limitations, but also at appropriate places in the manuscript (e.g. line 217ff, 376ff, 435ff).

The special section 4.2 “Study limitations” makes sense in this manuscript and addresses the issue (only briefly) but should not be disconnected from the results and discussion. I highly recommend referring to section 4.2 throughout the text and/or change the strong conclusions (like L.369 “Therefore, it must be the shallowness of the fjord at this point (indicated by low S) which governs the dynamics here.”) to more weaker expressions and elaborate a bit more on line 435ff.

Good point; we will make sure that the assumptions and limitations are more clearly articulated throughout the manuscript. Note, however, that the quote above is placed in context of comparing results of our artificial fjord settings where dS is the dominating control on retreat against settings where S is the governing factor. With this specific phrasing we therefore do not intend to draw any general conclusions about real-world glaciers, which would require a more detailed discussion on limitations, at this specific point in the text. We will clarify this distinction in a revised manuscript. Regarding line 435ff, we suggest to add a short explanation that not resolving vertical shear is less accurate where the bed slope is relatively steep, but that we do not expect our results to be hugely influenced as is pointed out in the comparison with a full-Stokes simulation above.

I assume you don’t see any advance during your simulations since you set-up the experiment to grounding line retreat. However, I suggest including a comment (in section 4), how and if your finding could also capture/predict GL advance? If your relationship of GL retreat and wetted area (or velocity) is to be applied in real-life, those areas might also have seen GL advance.

These are good suggestions that we are happy to address in the revised manuscript. Indeed, we do not see any grounding line advance, but previous studies ( Åkesson et al., 2018; Brinkerhoff et al., 2017) have shown that fjord geometry induces hysteresis in the retreat-advance cycle, where readvanceis inhibited if a fjord is widening or deepening ahead of a glacier. We expect this to hold for ourexperiments as well.

Specific comments

- L.6 Re-phrase, it is a bit confusing for an abstract: “We find that retreat in an upstream widening or deepening fjord does not necessarily promote retreat, but conversely, […]”

Rephrased as: "We find that retreat in an upstream widening or deepening fjord does not necessarily promote retreat, as suggested by previous studies. Conversely, ..."

- L.65 “larger suit of experiments” – why not mention you do 21 simulations?

Rephrase as: "Here, we use a numerical ice-flow model resolving two horizontal dimensions, we include a suite of 21 experiments and present a systematic approach..."

- L.143-154: The table 2 is very clear, but the text is partly confusing. For example: “we test 20 fjord geometries […] Additionally, we test […]” Clarify this paragraph in term of total numbers and subsets.

changed "additionally" to: "In the eight simulations outside our core experiments, we test asymmetric and longer perturbations to verify..."

- L.206 Define “very slowly” and “retreats quickly”

We define very slowly as <100 m yr⁻¹and quickly as >500 m yr⁻¹. This will be added in the next version of the manuscript.

- L.226: Please re-phrase for clarity: “However, in fjords that have a smaller S at xC than the reference fjord (bottlenecks and bumps), stable positions are also found where the fjord is narrow or shallow (small S). Therefore, S is also an important control on GL retreat.”

Rephrase as: However, glaciers in narrower or shallower fjords than the reference fjord (bottlenecks and bumps) can also temporarily stabilize where S is small. This shows that the wetted area constitutes an additional important control on GL retreat.

- Fig. 8a: Maybe adjust present day glacier front line colour…

Done.

- Fig. A1 e&f: include line for bed topography (like in Fig. 3 c,d)

Done.

References:

Akesson H, Nisancioglu KH, Nick FM. 2018. Impact of Fjord Geometry on Grounding Line Stability. Frontiers in Earth Science 6. doi:10.3389/feart.2018.00071.

Brinkerhoff D, Truffer M, Aschwanden A. 2017. Sediment transport drives tidewater glacier periodicity. Nature Communications 8:1–8. doi:10.1038/s41467-017-00095-5.

Favier L, Gagliardini O, Durand G, Zwinger T. 2012. A three-dimensional full Stokes model of the grounding line dynamics: effect of a pinning point beneath the ice shelf. The Cryosphere :12.

Goldberg D, Holland DM, Schoof C. 2009. Grounding line movement and ice shelf buttressing in marine ice sheets. Journal of Geophysical Research: Earth Surface 114. doi:10.1029/2008JF001227.

We thank the reviewer for the thorough assessment of our manuscript! We appreciate the constructive comments as they highlight where we need to justify our modelling choices more rigorously. We will address the comments as outlined below.

Methods

Why include an elevation-dependence for the flux at the inflow boundary and why have accumulation only in a small area upstream rather than a constant accumulation rate on the entire geometry as is typical in idealised modelling studies? Including the elevation dependence adds an unnecessary complexity that does not provide any additional insights.

The elevation dependence at the influx boundary is a way to parameterize the surface mass balance - altitude feedback. We acknowledge that including this feedback is not strictly necessary for the idealized experiments conducted here. However, we believe there are some good reasons for doing so: First, this feedback is well-established in glaciological theory  (e.g. Harrison et al., 2001) and has been shown to play a vital role in glacier and ice-sheet evolution (e.g.  Åkesson et al., 2017; Boers and Rypdal, 2021). We believe that not including it is a clear misrepresentation of reality and that the realism gained by including it outweigh the slight loss of simplicity. It is correct that we have simplified many of the complex processes governing glacier dynamics to reduce the number of degrees of freedom for the interpretation, but this is because many of these processes are poorly understood, badly resolved in models or costly to simulate. For the surface mass balance - altitude feedback, this is not the case. Second, the variation in influx is overall small compared to the total mass gains of the glacier. As is given by eq. 3 in the manuscript, the mass gain through surface accumulation amounts to 5.5 × 10⁹ m³ yr⁻¹, which can be compared to a mass gain of 1 × 10⁹ m³ yr⁻¹ through influx with a thickness of about 2 km at the influx boundary for our steady-state glacier. Towards the end of our simulation, the thickness at the influx boundary will have reduced depending on the fjord geometry, but not by more than ~1400 m (corresponding to an ice thickness of 600 m at the influx boundary). This minimum thickness translates to an influx of 0.3 × 10⁹ m³ yr⁻¹ which implies a reduction in total mass gains of only about 11 % over the entire simulation. So even under the assumption that including the parameterized surface mass balance - altitude feedback is unnecessary, it does not have a dominating influence on the retreat dynamics. Finally, it is worth pointing out that setting a constant accumulation rate on the entire geometry, as suggested by the reviewer, would also have led to a reduction in mass gain in the course of glacier retreat because the surface area of the glacier decreases as it recedes.

Why use the Budd-sliding law? This sliding law includes the dependence on bed elevation below sea level, which introduces additional complexity that can obscure the results (see also Brondex et al., 2017, for a study of the sensitivity of grounding line dynamics to the choice of the friction law). Where there are bedrock bumps or dips, this changes the basal resistance to flow.

The type of friction law is indeed one of the many model choices that inevitably will influence the results. Without an explicit subglacial hydrology model, basal drag needs to be parameterized. Many studies, including the current one, assume some dependence on bed elevation, either through the effective pressure (e.g. Morlighem et al., 2019;  Åkesson et al., 2021)  or through an elevation-dependent basal friction parameter (e.g. Aschwanden et al., 2019). This is true for Budd-type laws (where the basal effective pressure is given as N=ρ_igH−ρ_wg max(0,−z_B) ; cf. eq. 2), as well as for some Coulomb-type laws.  There is some physical rationale behind this elevation-dependency, whether implied through the effective pressure (as in the current study) and/or subsumed into a friction parameter (e.g. Aschwanden et al., 2019; Åkesson et al., 2021); ice generally flows faster further downstream (low elevations) than upstream (high elevations), and potentially weak sediments are more likely to be present in low-elevation areas. Again this will be a trade-off between real-world aptness and idealized simplicity.

The Budd law is one of the most widely used in the literature and we therefore think that our findings will be relevant to the wider glaciological community. A comparison across different friction laws is beyond the scope of the current study, and has indeed been done before in a slightly different setting, as pointed out by the reviewer. Nevertheless, we will follow the reviewer and include a more thorough discussion on the choice of friction law and the potential biases introduced. Specifically, we will clarify that the Budd law introduces an elevation-dependency, through the parameterization of the effective pressure, that is particularly critical for the fjords with bumps.

Use of different forcings to trigger retreat: again, this makes it difficult to compare different results. Why not use one (strong) melt forcing for all cases?

As is described in Section 3.2, we see the magnitude of forcing to trigger full retreat as a source of information, rather than a hinder to compare the results. This is because the forcing strength is one of two metrics to compare how efficient the different fjord geometries are in influencing glacier retreat (the second one being the residence time of the grounding line in a position of intermittent stability). For instance, had we simply chosen one strong melt forcing, we would most likely not be able to determine which of the small depression, medium depression or large embayment provide most stability, as the glaciers in these fjords all retreat after about the same residence time (Fig. 4). Furthermore, to choose the same strong melt forcing for all geometries would probably have induced faster retreat dynamics in those fjords where the glaciers also retreat with a smaller forcing. This would have reduced the level of detail and hampered our insight into the retreat dynamics. Finally, real tidewater glaciers react quite promptly to changes in ocean forcing (Khazendar et al., 2019). If, for instance, a slightly strengthened ocean forcing triggers glacier retreat, this will unfold immediately and not be delayed as the forcing strengthens further. For a study like ours, where we want to mirror this behaviour, this implies that it is the most realistic approach to assess glacier response to a forcing that is as small as possible. Therefore, we have chosen to try multiples of the reference forcing to trigger full retreat, and go with the smallest possible value. In our view, this is more realistic than setting one very strong value for all glaciers, because a real glacier would have reacted to a warmer ocean long before that ocean has warmed by a very large amount.

Lines 160-163 mentioning that doubling the melt rate leads to a reduction in calving and the grounding line is mostly stable. Again, such a choice of calving law is unfortunate, as it is not clear which dynamics are due to the topographic controls and which are due to feedbacks between melting and calving (see also Schoof et al.; 2017 and Haseloff & Sergienko; 2018 for discussions about how the choice of the calving law can alter grounding line dynamics).

We acknowledge that the choice of the calving law is an important control on the grounding line dynamics and calving front behavior. As discussed in the studies mentioned above (Schoof et al., 2017; Haseloff and Sergienko, 2018) , a calving law using a prescribed ice front position or a prescribed ice shelf length may have produced different results. However, both of these rather idealized approaches seem unsuitable as we have no such information for the future evolution of outlet glaciers. In fact, a comparison between different calving laws for Greenland outlet glaciers has been done before (Choi et al., 2018). One of the laws tested was a height-above-buoyancy criterion which equals the calving-at-flotation criterion discussed in Schoof et al. (2017)  if the tuning parameter q is set to 0. However, it was found that this law is not equally well able to reproduce observed retreat patterns as the von-Mises calving chosen in this study. We thus rely on those results in that we choose the calving law that was found most appropriate. In the absence of a universal calving law and considering the scope of our study, we find this approach to be the best available option, even though we are aware of the limitation that comes with it. Note that the calving stress threshold (σ_max in eq. 1 in the manuscript) has been set in accordance with typical values for Greenland outlet glaciers, as is explained in the manuscript (line 95). In the revised manuscript, we will more clearly stress the dependence of our results on the chosen calving law. 

Use of wetted area: I’ve never come across this term before and have difficulties subscribing to its usefulness (its application in figures 6 & 7 is highly doubtful, see below). What is wrong with just using the amplitude of the perturbation?

Indeed, to our knowledge, this term is not common-place in the glaciological context, although it is used by Catania et al. (2018) as `submarine area' (see their Fig. 4c). If it is deemed beneficiary by the reviewer or the editor, we are willing to change the wording from 'wetted area' to 'submarine area'.

Our wording was adopted from civil engineering and hydrology, where the wetted perimeter and the wetted area are commonly used to describe the size of the intersection between water and another media, for example the contact area between water and the hull of a ship (De Marco et al., 2017). Generally, the main advantage of using such a parameter (whether called wetted area or submarine area) as a metric for fjord geometry is that it integrates information about both depth and width of a fjord into one parameter. Our aim here is to find a universal relationship between fjord geometry and glacier response to external forcing; to that end, it is unavoidable to define such a quantity (we will answer below why we do think that its application in Fig. 6 and 7 is meaningful). Even when we designed the fjord geometries themselves, we explicitly wanted to introduce geometric perturbations of similar magnitude, but different type. Without the wetted area it would not be possible to define what a similar magnitude is when comparing basal perturbations with lateral perturbations. Only with the wetted area, there is a transparent quantitative link between, for example, a medium-sized bump and a medium-sized bottleneck. For real-world glaciers, the wetted area is also useful because it is not always possible to distinguish between basal and lateral perturbation in fjords with complex topography. Meanwhile, the wetted area can always be measured, provided sufficient information on fjord bathymetry is available.

Compared to taking the entire cross-sectional area (a `flux gate') of the fjord - glacier contact (the wetted area plus the area where the ice is above the water line), the wetted area has the advantage that it does not require information on frontal thickness. This implies that the wetted area can be calculated at any point in time and space for any fjord, regardless of whether a glacier is present at the moment or not. Clearly this shows the convenience of using the wetted area for paleo, present and future studies alike.

Results

Line 203 onwards: The terms “stable” and “unstable” refer to steady states, not transients, and are incorrect in this context as steady states are only attained at the beginning and the end of the simulation (see e.g., Strogatz, 2018). This needs to be rewritten to use appropriate terminology.

We follow the reviewer here and propose to use the term ephemeral grounding line position for what was previously called unstable grounding line position, and stagnant grounding line position for what was before referred to as stable grounding line position. This will be changed in the revised manuscript.

The presentation of retreat-results in figure 3 is not ideal as it is difficult to identify the important information from this plot. Better plot grounding line position in the center of the geometry vs. time and include both the results for the reference plot with forcing and the small, medium, and large perturbations. This should make it much easier to see where the retreat is fast and where it is slow and how patterns change with different topographic perturbations.

If we interpret the comment correctly, the reviewer would like to see a plot like Fig. S1, instead of Fig. 3 in the manuscript. Indeed, this has the advantage that several magnitudes of one perturbation type can be plotted in one subfigure. However, we prefer Fig. 3 as it is, for the following reasons: First, it demonstrates to the reader that we used a 2D model, as opposed to previous studies using 1D flow-line models only ( Åkesson et al., 2018). This underscores one of the novelties of our study. Second, we think that the current Fig. 3 allows a more intuitive interpretation of the retreat dynamics. At this point in the study, the reader has only seen a sketch of our study design (Fig. 1), and a visualization of the variations of the wetted area and its derivative (Fig. 2). With the figure as proposed by the reviewer, the reader cannot easily draw spatial correlations between retreat dynamics and specific features in the fjord. For instance, we think that it would be difficult to understand that grounding line retreat in embayments slows down where the fjord narrows in the downstream direction, because it is not directly visible from Fig. S1 at what location the fjord is narrowing. Therefore, it would be harder for the reader to gain an intuitive understanding of our results. Finally, Fig. 3 in the manuscript has the advantage that it shows processes which are very important for the retreat dynamics, and which can not be depicted by just plotting the grounding line position against time. For example, this refers to the ungrounding in the central part of depressions, a vital part of how retreat is `revived' after a period of grounding line still-stand. In summary, we do see the benefit of plotting a figure as suggested by the reviewer, but not at the expense of the current Fig. 3. Instead we suggest to add this as a supplemental figure in the Appendix. Furthermore, we will clarify in the caption of Fig. 3 that the spacing between the lines shows the retreat velocity, with lines that are closer together showing a slower retreat.

Some of the transient results are interesting and maybe counter-intuitive, but the presentation of the results and the unnecessarily complicated model assumptions make it difficult to trust those to be robust.

We are aware that our results may, to a certain degree, be influenced by the choice of modelling parameters, as is every modelling study. Besides choosing representative values where needed and applying commonly used parameterizations, we therefore dedicated Section 4.2 to how these choices may influence our results. This section will be extended in a revised manuscript, and references to this section will be inserted at appropriate locations in the text, following the comments by both reviewers.

Figure 6 & 7 and related discussion: Isn’t this simply showing mass conservation? For lateral variations, plotting Q/S is a proxy for the width-averaged velocity, which must increase where the geometry narrows simply due to mass conservation arguments.

Mass conservation will indeed play a role in the retreat dynamics, but it is far from the full story. First, if Fig. 6 in the manuscript was showing mass conservation, the highest (lowest) width-averaged velocity would occur where S is minimized (maximised), not dS. Crucially, using dS rather than S as the predictive variable (x-axis in Fig. 6), clearly supports our interpretation that the along-flow change in fjord geometry, not only the absolute depth or width, is an important control on grounding line retreat. Therefore, we believe that Fig. 6 is not simply showing mass conservation. To further cement this claim, we plotted v_GL over S (Fig. S2). Mass conservation would predict a strong relationship with a high (low) v_GL where S is small (large). This is because a glacier needs to speed up (slow down) when the fjord is narrow (wide) to maintain the same ice flux. Figure S2 shows that such a relationship does not exist. In fact, there is a weak tendency towards smaller v_GL for higher S, but this is not enough to explain the retreat dynamics that we observe.

Second, Q/S is not necessarily the width-averaged velocity, as is described in line 306ff: "Also, note that the GL flux is the product of the velocity v_GL and the flux gate area at the GL A_GL, that is Q_GL = v_GL x A_GL. The ratio Q_GL/S is thus proportional to v_GL when there is hydrostatic equilibrium at the GL (because in that case, S=0.9 x A_GL), ...". When the glacier front is grounded, there is not necessarily hydrostatic equilibrium at the GL, and then Q_GL/S is not directly linked to v_GL.

Finally, Fig. 7 shows the grounding line retreat rate dGL plotted over the wetted area S. It is not obvious to us why that would be related to mass conservation in a straightforward way.

As correctly stated in equation (3), with this choice the integrated accumulation at the grounding line depends on the ice thickness at the inflow and the prescribed parameters only, i.e., is constant over most of the domain. The width-averaged ice flux at the grounding line at the beginning and the end of the transient (when presumably a steady state is attained) should thus only differ due to differences in ice thickness at the inflow boundary. For transient model results the picture is less clear, but the dynamically interesting quantity is the width-integrated grounding line flux. Does that show deviations from expected steady-state results (ideally in simulations without the elevation-accumulation feedback)?

We are not entirely sure what the reviewer refers to here. The grounding line flux does indeed show deviations from steady-state values; in fact, this is a precondition for any retreat of any glacier where the surface mass balance does not change significantly. It is not clear to us why the grounding line flux would be expected to not do so, as is implied in the question. Unfortunately, we can not present any results without the elevation-accumulation feedback. However, we can show that the grounding line flux Q_GL varies considerably as the grounding line retreats (see Fig. S3 below). These variations are much larger than the decrease in influx at the upstream domain boundary (max. ~0.7 × 10⁹ m³ yr⁻¹, corresponding to 11 % of total mass gains, c.f. answer to first comment). Hence it should be clear that the large fluctuations in grounding line flux occur due to the variable fjord topography. 

Figure 8 & discussion: need to plot Q/S against dS to show that relationship still holds.

As is stated in the manuscript (line 329): "Plotting all available data points for Q_GL/S over dS at Jakobshavn, we do not find the aforementioned geometric relationship." In the original manuscript, we do not claim that the relationship Q_GL/S over dS (as shown in Fig. 6 for our idealized fjords) holds in a quantitative way for Jakobshavn. Therefore, we do not see the need include a plot of Q/S against dS. We do mention, however, that the relationship holds in a qualitative way, "such that an increase in dS is generally associated with a decrease in Q_GL/S and vice versa" (line 338). This `qualitative' support is shown in Fig. 8e,d. In response to the reviewer, we will rephrase this paragraph to make it clearer in a revised manuscript.

References