General comments

In this study (Part 2 of two related manuscripts by C. Schoof), a sophisticated mathematical model is formulated and used to compute the transient behaviour of the cavitation process and cavity evolution at the glacier bed. The model assumes a viscoelastic ice rheology and sets up boundary conditions in a way that allows the water thickness (h) to evolve seamlessly in space and time across contact transitions at which the ice lifts from the bed or vice versa. The approach can capture spatial variations in basal effective pressure (over distances within what the author calls the "process scale") and simplified hydraulics of water transfer between cavities and to newly cavitated areas. While the equations are formulated in three dimensions, numerical solutions are computed for a two-dimensional system, and, as in Part 1, Schoof prescribes the effective pressure at a fixed position along the bed to mimic its direct link to the ambient drainage system. The system is used to explore various behaviour in the dynamics of basal hydrological connections and disconnections.

A key motivation for the study – and a focus of its analysis – is to query some assumptions or parameterisations conceived or used in large-scale subglacial drainage models to describe evolving connectivity. Schoof’s simulations produce insights on this matter by showing a range of transient behaviour that exposes the limitation of the existing recipes. They also give a concrete demonstration of instances when borehole pressure measurements cannot adequately or usefully sample subglacial drainage conditions, highlighting the problems of interpreting such measurements. While Part 2 reveals new behaviour that occurs only in the dynamic system, some results are related back to those in Part 1 to show consistency between their modelling approaches.

The work in Part 2 clearly advances upon past modelling efforts and shows considerable prowess in both mathematical and numerical methods. It will be a useful reference for studies trying to improve the physical formulation of large-scale drainage models. The results should interest glacier hydrologists, both modellers and field scientists, so the subject of the manuscript very much suits the audience of The Cryosphere.

As for the Part-1 manuscript, I applaud the design of the numerical experiments in Part 2 and the methodical description of their results. However, while this manuscript is very readable, generally I find its text to be less refined than Part 1, with many glitches and typos, and a loose end in the physical explanation (see Specific Comment 1 below). Some passages are also imprecise (see my minor comments). I think that suitable revision to address these is necessary before the manuscript is accepted for publication.

Specific Comments

1. Oscillations in h (Sec. 3.2 onward; Figs. 3-5): you describe the “overshoot” and “undershoot” oscillations and analyse the observed factors behind their amplitude and decay rate with a good level of detail. Also you refer to the role of changing bed slope of contact areas (p. 15) and later discuss the implications of the oscillations (p. 24). However, I think that the physical cause of these oscillations is really never made clear or properly discussed in this manuscript. The sentence on p. 15 “These variations in normal velocity are presumably the reason for the significant oscillations in h-bar” doesn’t satisfactorily address the cause. Can you please fill this gap by adding a passage or paragraph --- at least to discuss candidate mechanisms if the correct one is difficult to determine? (Probably what looks like ‘propagating waves’ on the ice-base topography in Fig. 7 can aid the discussion.)

2. In Section 3.2, which presents the highly interesting "dynamic" run results in Fig. 3, it would help readers if you add a Supplementary Movie to accompany the figure and its textual analyses, such as those on p. 15. Since you made Fig. 2, going further to make a movie shouldn't be much more difficult. I leave this choice to you but I think that a movie will embellish the study.

Technical corrections: typos, minor suggestions, etc.

p1, line 13, "pressureized"

p1, line 17, "possibly other variables that can be computed by a large-scale model". This is vague. At least give an example.

p1, line 21, change “an average” to “a spatial average”?  I think this helps contextualise your subject

p1, line 23-24: the context is clearer if you insert the phrase "in the friction law " in the sentence "By contrast, basal water pressure is generally not assumed to be heterogeneous."

p2, lines 7-8. "pr"?  Suggestion: “The model *of* Rada and Schoof”

Eqn (2): correct the punctuation

p2, line 24-26: "… study instead how cavities can expand dynamically along the ice-bed interface from an access point where water is injected at prescribed pressure by an ambient drainage system". Clarify whether you're thinking in two or three dimensions. The next sentence specifies the number of dimensions, but that doesn't help us picture the idea of the current sentence.

p2, line 29: "varies slowly enough *in time*" -- this addition would make it clearer

p2, lines 30-34: your recount of the key findings of Part 1 here comes across as rather imprecise or vague, e.g.,

    line 30 “If cavity enlargement has occurred previously and cavity size has shrunk subsequently”. I can imagine that a cavity on a connected lee side that grows slightly and shrinks slightly, without extending over a bump top, also falls within this description.

    line 34 “reconnecting to an existing cavity is easier than creating a new cavity”. You probably mean a particularly kind of new cavity, not a new cavity that grows on a connected lee side as N decreases to below some *high* threshold value (e.g. N* = 8 in Part 1).

Fig. 1: (i) improve size of the arrow for h_w and the placement of h_w; (ii) in the caption, you should add a third sentence to say something along the line of "In this figure, the large cavity meets/overlaps with the stretch P, so it is connected to ambient drainage and its effective pressure is equal to ... [and so on]".

p4, nu is used here for Poisson's ratio but also later (p8 onward) for the small parameter in the shallow approximation

p4, around Eqn (3): I think that adding one or more suitable reference for this rheology (chosen for the ice) is necessary

Sec 2: To assure readers that the choices of rheology are sensible for the physical problem, I suggest that somewhere in this section you briefly explain why water compressibility (bulk elastic modulus about 2 GPa) can be ignored, while a compressible rheology is assumed for ice (bulk elastic modulus of 8-9 GPa; e.g. Table 1 of Neumeier (2018)), despite the stress coupling across ice--water interfaces. The reason probably is trivial and involves the very different dimensionless Maxwell times of the materials (i.e. when accounting for viscosities), but there may be other reasons.

p5, line 9, lower boundary *of the ice* (useful clarification, since b + h locates the upper of the two interfaces in Fig. 1)

p5, line 20-21: "Normal stress... , as water forces its way...". I suggest rewording this sentence because it isn't clear whether the "as"-phrase presents a scenario or reason.

p6, line 1: "impermeable except in specific locations at which water from an ambient drainage system can enter or exit the ice-bed gap". It would be useful if you describe explicitly (give actual examples of) what such entry/exit routes entail in this three-dimensional formulation. It is hard to picture a connection without knowing which direction or what materials are involved.

p6, lines 2-3, "for the remainder" isn't clear and you should "outside P" is that meaning is intended

p6, line 17 and Eqn (11d): is the correct symbol k or kappa?

p7, line 29, physics *of* (?) ice-bed contact areas

p7 (Eqns 11 e, f & g & Eqn 12): all sigma's and p's in this formulation differ from those in Part 1 where they had cryostatic overburden subtracted. I think that you should point this out in this section (even if any of the later analysis employs the subtracted version).

p8, the equations on this page lack numbering. Is this deliberate? Please check the journal's formatting guidelines.

p8, in the final scaling relation, it may be better to symbolise the water thickness scale by [h_w], as h symbolises the interfacial elevation (which is treated in the second-last scale relation).

p8, line 13: "defined" (towards end of line)

p8, line 14: by "forcing", do you mean "ambient"? Consider writing “ambient water water (which is used in this study as a forcing factor).

p9, Eqn 13a and preceding line: as mentioned for page 4, here you seem to be using nu for both Poisson's ratio and the small "shallow" parameter

p9, Eqn 13f: my attempt to derive this gives u-bar and v-bar instead of u1 and v1 in front of the derivatives. Please check.

p9, is there an Eqn 13g?

p9, Eqns 15c and 15d and next line: the conditions here seem to switch back into dimensional terms (for p_w at least), which comes across as confusing; that is, the p_w here doesn't seem to be the p_w in (13e), which I think is dimensionless. Please check.

p9, Eqn (13e) for sigma_33 at x3 = 0 seems to conflict/overlap with Eqn (14) (applied also at x3 = 0). Perhaps (13e) is replaced by (14) and/or it doesn't apply everywhere along x3 = 0?

p9, line 24, tau_b* --- you wrote earlier that asterisks are dropped

p10, awkward on lines 22 and 32 where the text switches back to referring to dimensional quantities when describing the numerical method of solving the dimensionless model of the last page

p11, line 14: the description here "code is implemented for both two- and three-dimensional domains" is a little confusing as the next line indicates that the code isn't used for three dimensions. The difference between "implementation" and "use" isn't clear.

p11 line 22: on declaring these choices for a and h0, it is useful to say that they make the N (dimensionless) in this manuscript directly comparable to N* in the Part 1 manuscript, as the effective pressure scalings are then the same. Section 3.1 later doesn't clarify this matter when comparing Part 1 and Part 2 results.

p11, line 23: hiccup after "In that"

p12, line 2: if I have guessed the intended sense here correctly, I would expect to read "highest" rather than "lowest" in this phrase. Please check.

p12 line 12, spurious curly bracket  [note: I’m counting downward from line 5]

p12, line 13: this lead phrase (“As measures… that … “)doesn't seem grammatical   [ Again counting downward from line 5 ]

p13, lines 3-4, while I understand this opening sentence, it would help readers if you add a sentence or insert a phrase to clarify whether Fig. 2 shows solutions in the moving or absolute frame of reference

p13, line 11, unclear what "the latter" refers to; clarify

Fig 2 caption, line 2, the phrase "the bed b is shown in grey" confuses b (the bed surface) with the bed interior (described as grey in colour)

Fig 3 caption, line 1 "effetive"

p15 & Fig 3:  perhaps this will be said later, or I've missed it. Although your focus on p15 is on the oscillations, it is useful to point out that the asymmetric response in Fig. 3 (h-bar doesn’t stabilise towards the same final value when N is step-changed to a certain value from different directions in this run) is related to the "irreversibility" of new cavity formation reported in Part 1 for the partially permeable case. This is in contrast to the reversible behaviour in Fig. 5 (fully permeable).

p15 last paragraph: you caution about the nature of the simulated oscillations at lowest N. But elsewhere in this section, you don't explicitly say whether you interpret the simulated oscillations at higher N (in fig. 3 and later figures) to be 'real', not dominated by numerical artifact --- although the writing seems to imply 'real'. Please clarity as a suitable place.

p16, lines 10-12 (irrelevance of viscoelasticity in Fig. 5): I have been wondering about this when reading p13-15. Can you please clarify whether viscoelasticity is also insignificant in the runs in Figs. 3 and 4 (besides 5) --- in causing the oscillations --- if that is true? 

Fig 4b panel: to help readers, please add the labels "cavity" and "contact", as done in fig 3; one or two places would do

Fig 5b panel: to help readers, please add the labels "cavity" and "contact", as done in fig 3

Fig 6 caption, line 2: is "18" a typo?

p19, lines 5-7 (delayed/final rapid increase in h-bar): Unlike the earlier phases of the evolution, for this final phase/part of h-bar rising, you don't give or hint at any physical mechanism. What controls or causes it? Or what delayed it, causing it to lag behind the rapid rise in theta-bar in Fig. 6a? Does the cause involve water transfer?

p19, line 14: in this passage it is worth pointing out also the brief recontact seen in panels g and h

Fig 8d: most steps in N have vertical lines. Add vertical line for the step at t = 260?

Fig 8 caption, line 2: hiccup in P value. Last line: I suggest moving "at t = 260" to elsewhere in the sentence

p20, line 1: columns of 9? figure 9?

Paragraph across p20-21: this description seems brief for the interesting result in column 1 of figure 9. If I'm reading Fig. 9a correctly, the connected cavity is longer (larger?) when water pressure (effective pressure N) is lower (higher)? Is this phase relation due to a time delay originating from viscous flow? Can you venture to say more?

p23, line 6, the value here (1.0653) differs from that in Fig. 2a-b

p23, line 10, "will also";   "will" seems redundant

p23, line 14-15, the message delivered here is "the subsequent growth of mean cavity depth h-bar... and of the cavitation ratio theta-bar ... causes a hydraulic connection to be established", but I don't think that it makes physical sense to consider these as cause and effect. (The next sentence seems to be fine as it uses the word "predictor", which conveys a correlation, not physical causation.)

p23, 2nd and 3rd paragraphs: these paragraphs seem to be written to address the context that (/the question whether) a specific variable threshold can be used (in macroscopic drainage models) as proxy for connection. These paragraphs will work better if you outline or remind us of the context at their start; doing this will serve to help the whole section. Currently this context emerges slowly, and I have long forgotten it since Sec. 1.

p23, line 25: "having a simple critical value hc" -- for what purpose?

p23-24: on these pages, you should highlight that here you're attempting to derive insights for drainage modelling in (I presume) three-dimensions from simulated behaviour in two dimensions. I am not sure that this translation from one to the other necessarily applies; the text on these pages conveys it as automatically valid for all aspects being considered. (This issue is linked to – but not the same as – the general limitations of using a two-dimensional model.)

p25, line 9, here you refer to the shape and volume of an "isolated" borehole. Do shape and volune matter because we are considering a borehole that has closed at the top by ice deformation? Please clarify in the text

p26, lines 4-6: in this passage, what end-to-end connectivity means is obscure to me

p26, line 20, typo in "caviites"

Reference:

Neumeier, J. J.: Elastic constants, bulk modulus, and compressibility of H2O Ice Ih for the temperature range 50 K–273K, Journal of Physical Chemistry Reference Data, 47, 033101, https://doi.org/10.1063/1.5030640, 2018.

In this paper, the author presents a highly sophisticated model to describe the dynamic evolution of subglacial cavities over a bedrock which can be impermeable over designated regions. The author carries out several numerical tests to show the ways in which cavity size and extent are related under different circumstances, and also explores the variation in effective pressure in isolated and connected cavities. These results are of great importance for the development of future large scale models of subglacial drainage. Therefore, I recommend the publication of this preprint after minor revisions.  

General comments

- I find the explanation of boundary conditions in section 2.1 after line 17 very difficult to understand. You say there are two alternatives for closing the system of equations for the viscoelastic fluid.  What I understand is that these two alternatives are either (10), which you enforce on P (or the whole bed when you consider a fully permeable bed), or (11e), which you enforce on the complement of P. In the case of (11e), the water pressure is given by the equations for the water column described previously. Is this correct? If so, I would consider rewriting this part of section 2.1 such that you first show the two alternatives (10) and (11e), and then explain how the water pressure p_w is modelled in the complement of P. I find this much clearer because what we need to close the equations for the ice dynamics are normal velocity and normal stress boundary conditions. 

I also have trouble understanding part of section 2.2, between lines 11 and 21. In line 14 you write: "imposes the boundary condition (14) only when (x_1, x_2) \in P is in a part of the bed to which the ambient drainage system has access". Does (x_1, x_2) \in P already imply that that point is on a permeable point and therefore has access to the ambient drainage system? I also see that condition (11e) is written as (13e) in the non-dimensional system. However, this is in conflict with condition (14). Shouldn't you include (13e) in (15)? 

- In equation (12) of page 8, you write far-field conditions for the basal shear stress. Previously, in Schoof2005, you enforced far-field conditions for the velocity. Why do include this new far field boundary condition here? I think it could be interesting to include an explanation for this choice of boundary condition in the paper. 

On this note, I am also confused about the sliding velocity variable u_b. You compute horizontal velocity perturbations u to this horizontal motion, yet you do not enforce the far field condition that u \to 0, right? In this case, the actual sliding velocity is u_b + u as x_3 \to \infty. I think it would be clarifying to mention this.

- In the numerical method, do you solve for velocity, stress, water pressure, cavity height and water height simultaneously? Or do you use any staggering of variables in time? I think it could help future researchers who wish to rexamine this problem to have access to the code you used. 

- The overshoot in the mean cavity size \overline{h} in e.g. Figure 3 is extensively commented throughout the paper. At some points you suggest it could be a numerical artifact, yet you refer to it to argue that cavitation ratio and ice-bed gap are not good proxies for each other (page 15, line 20). Therefore, it seems important to explore whether such oscillations are physical or numerical. This is obviously a difficult task and a full analysis of this phenomenon is out of the scope of this paper. However, a simple computational test would be to compute the cavity height after a step change in effective pressure for different meshes and time steps. If these oscillations were physical, we would expect the cavity height evolution to converge for decreasing mesh size and time steps. I suggest these computations be included in the paper, perhaps in an appendix. It would be very interesting to see this comparison for a case where dramatic oscillations occur, as in Figure 3 at t = 78. 

Specific comments

p1, line 13: "pressureized" > "pressurized"

p2, line 7: Unclear about meaning of "pr". Do you mean "i.e."?

p2, equation (2): Add full stop.

p3, line 10: "practical computational reasons" - What does this mean exactly? Excessive computational cost, unsuitable numerical model for 3D, difficulties in formulating/implementing computational tests? This should be clarified here. 

p4, equation (3): I think some readers might not be familiar with the mathematical description of an elastically compressible upper-convected Maxwell fluid. Adding a citation where equation (3) is derived/explained would be very helpful.

p4, line 13: close brackets.

p4, line 19: "then with the change in stress related to the corresponding linearized strain as (eq)". Rephrase this clause, it is phrased incorrectly. 

p4, line 28: Avoid initiating sentence with mathematical symbol.

p5, line 11: "ensures ensure" > "ensures"

p5, line 17: Consider rephrasing the sentence "First, the standard assumption in dynamic models of subglacial cavity formation (references) has been ...". Perhaps write "First, we consider the standard assumption in dynamic models of subglacial cavity formation (references), which has been ..."

p5, line 26: "in contact areas, normal velocity is prescribed". If by contact areas you mean areas where h = 0 (which is the most intuitive definition), this sentence is not correct. We will also have contact areas which are about to detach, \partial h / \partial t != 0. In this case we prescribe the normal stress and compute the normal velocity. 

p6, equation (11b): Add comma.

p6, line 21: "flux q" > "the flux q"

p6, line 21: Consider rewriting "... by the first, pressure-gradient-driven term". Perhaps "... by the first component of the flow, the pressure-gradient-driven term".

p 7, line 9: "implying that a source term that is omitted in (11a)" - this clause does not make sense, please correct.

p7, line 29: "also capture the physics ice-bed contact areas" - Typo?

p8, line 23: "N* is the (scaled)..." - Avoid starting sentence with mathematical symbol.

p9, equation (15b): Add comma.

p10, line 13: Specify that the  mixed finite element method is used to solve for the velocity and stress variables. A mixed FEM is used in Stubblefield2021 and deDiego2022 to solve for the velocity and pressure and the velocity, pressure and normal stress at the bed, respectively. 

p10, line 19: I can see how a moving frame eliminates the advection terms in (13a), since these are advected by (\overline{u},0). However, I do not see how the advection terms disappear in (13f). Do you mean (13b)?

p11, line 19: "transverse normal stress" > "the transverse normal stress"

p11, line 24: How small are the intervals around x_P?

p12, line 11: Indicate which endpoint is upstream and downstream for each cavity.

p13, line 1: "the inherent heterogeneity involved in an unstructured mesh" > In what way does the degree of uniformity of a mesh influence possible oscillations in the cavity shape?

p13, line 5: "at least for the moderate values of N for which the dynamic model produces a recognizable near-steady state within a reasonable time span" > Does this mean that for smaller values of N, the difference between the cavity shapes produced with both models start to differ visibly? If so, I suggest that an additional panel be added to figure 2 for a value of N for which the models in part 1 and 2 start to differ. If future work is to be produced on this topic, researchers should have an idea of the ways in each the numerical model you propose produces potential inconsistencies. Figure 2 as it stands now indicates an almost perfect consistency between both models, yet what you write suggests the contrary. 

p13, line 17: "steady state mean water depth \overline{h}" - \overline{h} refers to the mean cavity size, which coincides with the mean water depth in the cases considered here. I suggest you avoid refering to \overline{h} as the mean water depth here because it could be confusing for the reader. 

p14, Figure 3: It could be interesting to show values for e.g. \overline{h} at the steady states obtained with the model from part 1 if the same time history for N was followed (allowing for quasi-steady states to be achieved by small changes in N, as opposed to the step jumps we see in Figure 3). This would give a valuable insight into how the dynamic evolution of cavities differs from its steady counterpart, which is one of the main goals of this paper.

p14, Figure 3: This figure would be more readable if the ticks of the x axis were aligned with the vertical grid lines. Throughout the paper you refer to the times where jumps in the effective pressure take place (e.g. t = 78) and its not entirely obvious which points these are. The same goes for the remaining figures in this paper of a similar type. 

p15, line 24: "That contact area motions occurs around the top of the prominent bed protrusion at x = 0.8." - This sentence does not make sense, change "that" > "these"? 

p17, Figure 5: As I wrote above for Figure 5, it would be very nice to include results for the steady solution here too, in order to see for example whether the oscillations in \overline{h} occur around the steady states predicted in part 1.

p19, Figure 6: "bed 18" > "bed given by (18)"

p19, line 16: "the a less-advanced" > "a less-advanced"?

p19, line 1: "an initial advance of the cavity end point from c ≈ 4.1 to c ≈ 4.5 over a time interval around 10^-2" - I do not see this in figure 6. Over a time interval of 0.01 I see an advance from around 3.9 to around 4.1.

p19, line 4: "htis" - "this"

p19, line 13: "leading oscillatory" > "leading to oscillatory"

p19, line 22: "section part 1" > "part 1"

p20, Figure 7: This figure could be improved by adding visible marks indicating the endpoints of the cavities.

p20, line 1: "9" > "figure 9"

p20, line 6: "the smaller bed protrusion upstream of N" - Do you mean M?

p21, Figure 8, caption: "P = {./65}" > "P = {.65}"

p21, line 8: "a extended" > "an extended"

p22, line 9: "the greater ability of the solution to relax towards a steady state" - The reader could judge the validity of this statement if information on the steady states was included in Figure 9. I suggest values for N_M and x associated to steady states be included in Figure 9.

p22, line 10: "what the extent" > "what extent"

p23, line 17: "it is plausible a critical value h_c could plausibly be defined" > Avoid repetition of plausible/plausibly?

p 24, line 27: "Once a a set" > "Once a set"