Review of “SUHMO: an AMR SUbglacial Hydrology MOdel v1.0”

Submitted by Felden et al. to Geoscientific Model Development

General comments:

This manuscript documents the development of a new model that simulates the evolution of water flow and pressure beneath glaciers and ice sheets. SUHMO adopts the general equations and continuum approach of the SHAKTI model (Sommers et al., 2018) with modifications to facilitate incorporation of adaptive mesh refinement (AMR), in order to resolve individual subglacial channels.

Overall, the paper is clearly written and well supported with references and figures. I do not find the simulations based on SHMIP tests to be especially helpful as there is no AMR and channel resolution involved in the SUHMO results presented. However, I understand the motivation of the authors to include those results as a sort of benchmark or verification of the approach, and I leave the decision up to them whether or not to retain those.

In general, some aspects of the numerical experiments and results presented in the paper could be more clearly explained. For example, the melt rate in SUHMO assumes no geothermal flux and no frictional heat. These can be large sources of basal melt, particularly frictional heat for fast-moving glaciers, and the decision to neglect these should be explained.

As already commented by the Editor in the online discussion: the model code, input, and output need to be made available in a repository that complies with GMD’s guidelines.

Please see specific comments and technical corrections below.

With some revisions, I feel this manuscript should be published. As conveyed by the authors, challenges remain in incorporating the influence of subglacial drainage into ice dynamics models, and this AMR approach is a promising step to evaluate the role of individual drainage features and better understand the spatial resolution necessary to capture relevant subglacial pressure for large-scale ice sheet simulations.

Specific comments:

Line 11: Consider strengthening this last sentence by including the main point or points of what you find out about effective pressure as resolution increases (instead of simply saying you discuss it, which is vague).

Lines 72-73: Studies that have used subglacial hydrology modeling on real Greenland glacier domains find widespread areas of channelization (for example, Cook et al., 2020). In that light, I’m not sure that channels are only expected to occur in very localized areas – each channel is distinct, yes, but with AMR you would need to have many areas of refinement that could cover a large area of a glacier bed. You may want to consider re-wording this sentence.

Line 73: You could define AMR here, in the first instance of using it.

Lines 74-75: Suggested re-phrasing: “We follow the approach of Sommers et al. (2018), with adaptations to implement a subglacial hydrology model using the Chombo AMR framework (Adams et al., 2001-2021).”

Line 75: Consider briefly describing what Chombo is here.

Line 102: What value of omega is used?

Line 118: The term in the melt rate that accounts for changes in the pressure melting point with changes in water pressure as written here is not negligible in places with steep topography, for example (see Creyts and Clarke, 2010 about supercooling). There are good reasons to drop this term, however, that have to do underlying assumptions about the ice and water pressure being equal at the interface.

Equation (6): Why introduce the general b’ here and not in Eqn. (1)? I would recommend to either present a more thorough treatment of partially filled drainage, or just go with the assumption from the start that the gap being ice and bed is always filled with water.

Line 125: This method of using beta for the opening by sliding is based on GlaDS (Werder et al., 2013).

Line 130: Please describe the “creep length scale” and the physical justification for Equation (7). What qualitative behavior is captured with this equation? Why does creep shut off below a threshold? (Hint: it may have something to do with ice being supported by asperities on the bed). GlaDS and SHAKTI use the gap height (b) as this length scale. Why is that, and why should it be improved upon?

Line 131: Again, see comment above about b’.

Table 1: tau_b and G are both listed as 0. You should make clear in your results that the melt rate does not include contributions due to geothermal flux or frictional heat from sliding. Why?

Table 1: The description for omega is not very informative. You may want to include a brief description near Equation (2) about what this parameter is and how it functions to make the laminar-turbulent transition.

Table 1: Why did you choose this value of A that corresponds to very cold ice? The equations described assume that the ice and water are both at the pressure-melting-point temperature, so a larger flow law parameter would be more appropriate.

Line 151: The diffusion-like term is basically a diffusion of gap height. It may be worth explaining in more detail what this represents physically.

Equation (10): I would not consider the second term in the parentheses (the pressure-melting term) to be dissipation. I suggest naming it separately.

Eqn. (16): What is rho_c?

Line 211: If only one Picard iteration is used, what’s the point?

Equation (17) and Lines 214-216: I recommend using a different letter to indicate time in the superscript to avoid confusion with the flow law exponent n.

Line 231 and Eqn. (19): Here too, the use of little n as n_tot could be confusing.

Line 237: The domain (64x16 m) used for the channelizing test case is small. What about boundary effects?

Line 239 – 30 m3/s is pretty high flow for a moulin!

Line 240: Neumann (flux) boundary condition at the outlet and Dirichlet (value) boundary condition upstream? I think this is probably a typo and should be switched. Typically, the Dirichlet condition would be at the outlet to set head equal to bed elevation (for example, equivalent to atmospheric pressure). And the no-flux Neumann condition makes sense for the upstream condition based on your results.

Line 242: I was confused on first reading this, thinking that the Gaussian was referring to a time-varying input rate, whereas it was already stated that the input rate is 30 m3/s. It would be helpful to include wording that makes clear that the Gaussian distribution refers to the spatial footprint of the point source at the bed.

Line 248: It may be worth commenting here about why this type of super-narrow, tall opening should not be interpreted as a physically realistic basal crevasse.

Line 302: What do you mean by the effective system? Is this a typo for “efficient”, or some other intended meaning?

Lines 306-307: I thought the moulin input is specified as constant. It is not clear what is referred to as Gaussian here.

Figure 5b: How high does the discharge calculated by SHAKTI go in case A1? (The vertical axis does not extend sufficiently). This seems like an error somewhere (quite possibly in the SHAKTI results submitted to SHMIP).

Lines 337-343: 5180 m3/s into 63 moulins = 82 m3/s into each moulin. This seems very high and could benefit from better justification or explanation of why you choose to use extreme input. What happens with more realistic moulin inputs (not as high)? Do you produce channelization?

Figure 9 – suggest moving the x,y,z coordinate in panel a to the right lower side so that it is oriented where x=y=0 (the red dot locator helps orient, but I found it confusing nonetheless). “Total height” could be called “surface elevation” instead.

Figure 9 – panel b: is moulin input in m/s or m3/s?

Line 351: What specific criteria for gap height and melt rate were used to determine regridding?

Line 370: Great! This is a big deal and will help make the case for AMR. It would be helpful to include a number here of what the minimum resolution is that you find to resolve channelization.

Line 383: I am not sure that it is correct to say that melting of ice walls was discarded in the original derivation of the SHAKTI equations. The equations treat gap height as a one-dimensional quantity, so two-dimensional cross-sectional area of a channel was never part of that derivation.

Conclusions section: Make a point of quantifying the minimum mesh resolution that you find necessary for resolving channelized features somewhere in here as a main take-away point.

Lines 391-395: This paragraph is arguably true, but I suggest rewording to strengthen the case for resolving channels in SUHMO. Describe why this may be necessary and helpful in ice-sheet modeling. (The way it currently reads make it sound like coupling with BISICLES is motivated by making sure individual channels don’t matter – but we don’t know yet if that’s the case until you do it).

Technical corrections:

Line 3: comma before “however”

Line 71: typo in “subglacial”

Line 82: convergence analyses… are presented (or convergence analysis… is presented)

Line 124: should be kg m^-3 (missing a negative sign in the exponent)

Line 162: Insert word “accuracy” after algorithm

Line 173: missing space after “of”

Figure 1. Should the schematics on the bottom left and right have delta x₀ (instead of delta x) in the grid level labels for consistency?

Line 276: Period after Figs.

Line 293: Extra ) after 2a

Line 308: Extra )

Figure 9: The orientation of the two panels being different is spatially confusing. Is it possible to orient them so the red dot location is more consistent between perspectives? For example, rotate panel b 90 degrees counter-clockwise. Similarly for Figure 10b.

Line 366: Missing period after Fig.

Line 366: Either includ gap height b, or remove N

Line 373: Using lower-case n is again slightly confusing here.

The manuscript by Felden et al. describes the newly created subglacial drainage model SUHMO.  Its mathematical formulation is based on a linked-cavity inspired water-sheet drainage mechanism which has been used in many recent models.  Unconventionally, it allows for melt due to dissipation to occur and thus flow localisation is possible (akin to channelisation which I use interchangeably here), as was also done in the SHAKIT model (Sommers et al., 2020).  However, in an implementation such as in SHAKIT, such a flow localisation mechanism leads to unbounded localisation.  In a novel approach, SUHMO stabilises this localisation with a diffusion term, which renders the equations well posed.

The numerical implementation is based around the adaptive mesh-refinement (AMR) library Chombo; the ultimate aim is to eventually couple SUHMO to the well-know BISICLES ice sheet model which is based on the same library.  The manuscript presents the numerical algorithms, based on multigird strategies, used to solve the discretised PDEs in good depth.

The manuscript then presents comparisons to results of SHAKIT from the Subglacial Hydrology Model Intercomparison Project (SHMIP).  These show that the two models produce very similar results, except in a no-SHMIP test-case where channelisation occurs the two differ markedly.  Last the manuscript presents a test case in which its channelisation and AMR capabilities are displayed.

This is a well written manuscript which is well suited to be published in GMD once a few changes have been implemented.

Please note that I am not an expert in AMR nor multigrid methods, thus I will not comment on those parts of the manuscript.

Major points

Note that I have quite a few major points listed here.  This is merely a reflection of my interest in the model and should not be seen as a harsh critique of the manuscript.  In fact the manuscript is well suited to be published without implementing these suggestions.  And indeed some of the suggestions may go beyond what should be expected from a GMD paper.

*Show-case channelisation capabilities*

Whilst the comparison to the SHMIP results of the SHAKIT model show that SUHMO indeed can behave very similarly to that model (as is expected), it would be more informative to show where there are differences.  Thus I would prefer to see SHMIP experiments where there is significant channelisation.  Note that the SHMIP instructions suggest to tune to the GlaDS output of runs A3 and A5, with the latter being channelised up to mid-domain in GlaDS.  The SHAKIT submission only tuned to A3 and only shows channelisation in A6 (GlaDS has channels in A4 to A6), thus it is one of the models which shows the least channelisation in the SHMIP study.  Also, the shown suite-B is channelised in GlaDS but not in SUHMO & SHAKIT.  (Aside, presumably in a structured grid model such as SUHMO some noise would need to be introduced, e.g. via initial conditions, to nucleate channelisation in a setup like SHMIP-A.)  Also, showing a comparison to GlaDS outputs for SHMIP cases would be informative, as that is currently the canonical distributed+channelised drainage model.

In a similar vein, the purpose of SHMIP suite E was to look at effects of an overdeepening which are influenced by pressure-melting point effects, i.e. by setting $c_t>0$.  Thus it would be nice and interesting to see results of SUHMO results with pressure-melt term turned on.

All in all, I would ban the one-to-one SHAKIT-SUHMO comparisons to the appendix or supplement as it merely shows the correct implementation of the model and not its novel capabilities.  Replace these with results which show the channelisation capabilities of SUHMO in action.  I realise that this is probably a reasonably large undertaking, so I would not insist on it, but I do think it would make the impact of the manuscript much larger.  After all, one of the coolest bits of this model, in my opinion, is the regularisation of flow localisation with the novel diffusion term.

*Figures of channels*

Everyone and their dog want to see figures of channel systems.  On this the manuscript fails.  The closest gets Fig.11 which shows the channel in this setup but very small and blue-in-blue.  It would be cool to see the channel system clearly for this case but also for SHMIP A6 and E-Suite.

*AMR and channelisation*

The authors mention in Section 6 that for channelisation to occur resolution needs to be good enough.  It would be good to elaborate on this point some more.  Does the used refinement criteria allow to always capture the channels or could some be un-noted and thus not simulated?

Also, in the small-scale experiment (Fig. 4) the channel is resolved with ca. 10 nodes across its width.  How well resolved is the channel in Fig.11 for R4?  Should a channel always be resolved with several grid points across its width?  Or is one enough?  I think it would be ok to leave these questions open for a future study but mentioned they should be in the discussion.

*Regularisation and comparison to "proper" channels*

It would be interesting, and I think eventually necessary, to compare the SUHMO simulated channels to more classically simulated R-channels such as in GlaDS.  Do they conduct the same discharge at the same pressure gradient?  How do parameters translate between a classic R-channel formulation and this regularised sheet formulation?

Smaller comments

Please only use abbreviations if really needed, i.e. reducing the length considerably and reducing line-noise for the reader.  Thus remove the abbreviations "GrIS", "AIS", "GMSL", "GC".  I'm a bit unsure about the "AMR" in the title, albeit it should be ok for GMD.

One of the main aims of this model seems to be eventually coupled to the BISICLES ice sheet model.  I think this is worth mentioning in the abstract.

One of the main reason that subglacial hydro modelling is not further along is the lack of observations.  This should be mentioned in the Abstract and Introduction.

Many figures are too small to be legible.  For some of them, extending them to the full page width should be enough.

It would be good to look at convergence also for the case of section 6 as that features a channel.

I feel that the SHAKIT parameters used in SHMIP lead to too much sheet flow and not enough channel flow compared to the other models in SHMIP.  I don't know what the reason was to pick those parameters, maybe SAHKIT struggles with high channelisation or the parameters were just the optimal ones.  Irrespective of the reason, the SUHMO-results presented, based on the SHMIP-SHAKIT parameters, end up having to use very high inputs (both in run for Fig 4 and Fig 11) to get channels.

It would be good to have a master-scrip (or master make-file) which produces all the figures in this manuscript.  Maybe there is one already, but 2min of browsing the github repo did not turn up one.  State this in the "code availability" section.

Line by line

9: should it read "MultiGrid" to match abbreviation?

16, 17: citation needs an "e.g."

20: that "likely" is probably IPCC speak, why not reference the latest IPCC report?

29: also cite Gilbert et al 2022 https://doi.org/10.1029/2021GL097507

31: write "surface-lake water"

32-34: this sentence does not work, rewrite

37: compare length to lengths and not area!

39: cite Wertman (1962) for sheets, and Walder (1986) and Kamb (1987) for cavities

67: write "The SHAKIT model..."

67: isn't it "SHaKIT"?

75: write "We propose a small but significant ..."

84: "results"

95: I find "water gap height" odd as there is no gap in the water.  Maybe "gap height" or "water height" or "water thickness", "water-filled gap height"...

110: it would be nice to contrast this formulation to the more commonly used Manning or Darcy-Weisbach formulation.

119: ".. and dropped from some similar models..."

120: I find the "effective drainage-system capacity" a confusing term, something like "bed-ice gap" would be much clearer.  Irrespective of the term used, it should be clearly explained and its relation to b, the "water gap height".  (or as the other reviewer writes, just use the same variable without discussion)

Eq.6: this cavity-sheet formulation was first introduced by Hewitt 2011 https://doi.org/10.3189/002214311796405951, cite.

Eq 7: this needs an explanation.  Why do is this term introduced?  Presumably to avoid complete closure of the sheet?

153: "... dissipation term of the melt rate (Eq. 5)."

Eq 12: what are A, u and F?  Define.

173: missing space

Fig 1: illegible, should be full width.

202: it should be stated earlier that this system is solved.  Thus far it was only ever one equation.

Fig 2: illegible, should be full width.  Also note in the caption that (c) has a vastly different scale.

Fig 3: larger.  Why is P_w also plotted.  Is that not a simple dependence on h, and thus has essentially the same error?

242: sentence does not make sense, rewrite

272: in the spirit of SHMIP, A4-A6 should be channelising test cases.

296: these instructions are also in the supplement of de Fleurian et al., maybe better also cite that less ephemeral source.

311: The point of SHMIP-E is to see effects of overdeepenings and those effects show only with a pressure melting term.  Thus removing it is a bit odd, except to have a one-to-one comparison to SHAKIT.  See "Major points"

336: state what the smallest \Delta x is (24m).  Also, to resolve channels with several grid points, probably a R_8 would be neeeded (1.4m).  Maybe state that R_4 does not resolve channels?

Fig 9: far too small!

Fig 10: even smaller!!

370: write in this paragraph something on whether we can be sure that AMR will refine in the places it needs to, to capture all channels.

371: state that refers to an ice flow model

373: I don't think this notation has been explained.  Maybe better to just write it?

377: Iken & Binschadler 1986 (Fig 6) would be the canonical example,

no?

Fig 11: I feel this is the central figure of this paper and that it could go a long way to sell SUHMO well.  In this light, it needs to be better legible, and larger.  The channel in the top panels needs to be visible.  A scale bar would be nice.