This paper models the opening of hydraulic fractures below supraglacial lakes, and the subsequent spreading of the basal fracture, uplifting the overlying ice. The aim of the paper is to investigate the role that different rheological models for the ice have on the system. Results are compared with field data from Das et al. 2018 for a lake drainage event. The study finds that the inclusion of viscous creep in the ice is important to accurately capture something close to the observed behaviour: an elastic ice model does a bad job, even though the dynamics occur over a short timescale (a few hours) - presumably because of the large stresses involved.

This is an interesting study and I think the conclusions are useful for the community. However, there are quite a few points that I think could be explored in more detail, and the model results need a bit more exploration to be entirely convincing. I think it could use some fairly significant revisions, as outlined below, to increase its impact, and also to alleviate various issues of unclarity or inaccuracy.

Specific comments are below, along with various more technical comments.

Specific comments:

1) The problem is studied in two dimensions, and the authors go to some effort throughout to argue that this is a reasonable limit to consider. The other simple ‘end-member’ option would be a radial (i.e. axisymmetric) profile, as mentioned around line 75. I don’t fully follow the reasoning in the paper here: there are statements that don’t make sense, like “While our 2D model for the horizontal basal crack propagation and the basal uplift is valid for the axisymmetric assumption…” Presumably this should mean something like “The 2D planar model construction could be straightforwardly adapted to describe instead an axisymmetric spreading”. The point is, the construction is not complicated, but there are different- and non-trivial - geometric factors which change some details.

Fundamentally, it seems to me that it should be straightforward to do the whole problem in an axisymmetric geometry as well as the 2D planar geometry - because it uses the same ideas and is still mathematically 2D. And this would be very valuable, because ‘reality’ is somewhere between the two (planar and axisymmetric; although arguably it is closer to axisymmetric than planar) and so comparing solutions for each would greatly help the impact of this work.

2) The results in figure 7 are concerning as it stands. The curves - particularly the black curve - are curiously non-smooth, and the mechanisms / reasons for this are not at all clear. There is some discussion around line 388 about this, but the explanation is not very convincing, and much more evidence needs to be presented to convince the reader this is not some funny numerical artefact. I can’t see which aspect of the mathematical formulation is giving rise to this behaviour - for example, the black line goes flat for a period, and then increases suddenly. Is this robust to numerical resolution? What aspect of the model allows the crack to halt propagation for a period and then restart motion? What physics is controlling the length of time the crack is stationary for, and indeed, the time at which it decides to become stationary? Much more convincing analysis is needed of this behaviour.  The point is discussed again around line 450, but again I don’t see how the model is giving this ‘episodic’, almost stick-slip-like behaviour. Perhaps a plot of something like the pressure at the tip of the spreading crack would help to explain this phenomenon.

3) The comparison with data from Das et al. is interesting, and a bit more could be made of this. The main disagreement seem to be that the water-level change (i.e. the flux into the conduit) goes quite wrong in the model: much more water gets into the crack system than the model predicts. Interestingly, given the model is predicting the wrong amount of water in the system, the uplift prediction is quite good initially (although it also goes wrong at later times). The explanation about a pre-existing damage network seems plausible. I was surprised not to see more discussion about the possibility that the bedrock is not frozen to the ice: if the basal ‘crack’ or conduit can spread without cracking (a zero-fracture-toughness limit) then presumably the crack would spread further and allow more water in, without necessarily increasing the localised uplift (because the water has spread laterally further). It would, presumably, be straightforward to consider simulations with different fracture toughnesses for the ice-bedrock interface - and this seems valuable anyway, because we don’t really know what that value should be.

It also seems likely that the 2D planar assumption has quite significant errors as the spreading at the base continues, compared with an axisymmetric model, which is perhaps behind the later-time disagreement in the uplift (the geometric constraints are rather different for spreading as a circle compared to spreading as a line)? Again it would be useful to be able to compare the model predictions.

4) One of the aims of this work, as I understand it, is to highlight the role of the viscous rheology, and the fact that the interplay between viscous creep and elastic deformation can be very important in these processes. I think this point would be aided by a bit more analysis of the relative importance of the two modes of deformation. Specifically, the work compares the ‘linearly elastic’ model (just elasticity) with the ‘viscoplastic model’ (elasticity and viscous creep), but we don’t really learn about how important viscous creep is relative to elasticity in the latter. i.e., one might be tempted to conclude that the role of viscous creep is dominant here, and that a third model of pure viscous creep (no elastic deformation at all) would do fine. It would be interesting to look at how much of a role elasticity is playing in the model, to be able to draw a clearer conclusion about how much the interplay of elasticity and creep is important here.  That would be a helpful qualitative conclusion to draw from the work.

More technical comments

- Throughout: the non-linear visco-elastic formation is throughout referred to as a ‘viscoplastic’ law. I know this terminology is sometimes used to describe Glen’s flow law, but it isn’t strictly correct, and anyone from a non-Newtonian fluids / rheology background would be confused by its usage here. They would traditionally Glen’s flow law as a shear thinning viscous rheology - and the model used here would be a viscoelastic model: the ice comprises elastic deformation (recoverable) and non-linear viscous deformation (non-recoverable, and given by a shear-thinning model). A ‘viscoplastic’ model would typically be taken to indicate that there is a plastic ‘yield’ stress, below which there is no non-recoverable deformation; at which the material deforms plastically; and above which the material flows viscously (see e.g. much literature on visco-elasto-plastic models). I would favour not describing the formulation here as ‘viscoplastic’.

- p.6 Figure 3: needs to say this is the range of Maxwell times - the caption just says ‘range of time-scales’ which could mean anything.

- p.6 equation (3) has some weird typesetting on the third line (missing equals?)

- p.9 The authors choose a Manning-Strickler turbulent flow law following Tsai and Rice and others. Do any of the results have any appreciable dependence on this choice as opposed to other turbulent laws? (e.g. see opening of Dontsov 2016 J. Fluid Mech. ‘Tip region of a hydraulic fracture driven by laminar-to-turbulent fluid flow’ or final section of Hewitt et al. 2018 J. Fluid Mech. ‘The influence of a poroelastic till on rapid subglacial flooding and cavity formation’).  Those studies also give details of how to map to a laminar regime if the fluid velocity / crack dimensions become too small: could this be important at later times? (Particularly in the case where the linear-elastic model shuts off the water supply to the basal crack, which then slowly continues to spread - see point below.)

- p.10 Equation (12) does not make sense. The condition is supposed to be a fixed pressure (i.e. at the base of the lake), so what is this ‘penalty’ term added to (12), and how is the parameter chosen (the text says it is chosen ‘to be large enough to enforce this inflow condition’, which isn’t clear: it is part of the inflow condition, so the value you pick for it will change that condition.) In addition, if the problem is being driven by a fixed pressure condition at the start of the fracture, than I don’t see what the variable p is chosen to be at that location - it should surely be p = p_ext, but that seems inconsistent with equation (12). This means I don’t really know what is being imposed on the pressure in the model, which is awkward when trying to interpret the results of the model.

- p.11 equation (15) seems to have the wrong units for a heat flux (is there a missing c_p?)

- Equation (16-18) I don’t think eta is defined here. More importantly, the expression in (18) is obviously wrong and seems to have been lifted from elsewhere - the factors of pi and the power of (3/2) rather than (1/2) on the (t-t_0) indicate this is the corresponding expression for an axisymmetric fracture, not a plane fracture.  Equation (19) has the same issue - and so this error will propagate into all of the numerical results.

- p.15 Figure 6 looks a bit odd: it would be much more clear if the colours filled the gaps, rather than being a line in the centre of the gap.

- p.16 Line 370. The statement that the crevasse closure stops the propagation of the horizontal crack is inconsistent with figure 7, where the crack continues to spread slowly after the closure (i.e. the blue line keeps rising after about t=45min, when the gap closes). It is unclear to me why the lateral fracture continues to spread, even though the entry has closed over: if it has enough pressure to keep spreading at the tip, why does it not have enough to force the original opening back open? 

This is an interesting model study of an important process. The model represents an advance on previous treatments of hydrofracture, and has yielded some interesting new insights into the relative importance of thermal and deformational processes in controlling rates of hydrogfracture propagation and supraglacial lake drainage. The anonymous referee has raised a number of technical points regarding the paper, and I shall not repeat them here. With regard to terminology, I agree that the term 'visco-plastic' is not appropriate and should be changed to "viscous" or "non-linear viscous" throughout. It is also unnecessary to refer to the "so-called Glen's flow law" (line 126). Although some have used other names for this flow law from time to time, "Glen's Flow Law" is in very widespread use, so the "so-called" is superfluous.

My main comments concern model formulation and how it may relate to reality. The 2D plane geometry seems to me to be perfectly adequate to explore the issues of interest, and there is perhaps no need to add experiments with an axisymmetric geometry (full 3D would of course be much better still). More serious for present purposes is the assumption of a frozen bed. In Greenland, the bed is temperate below most of the ablation zone, and water is certainly present at the ice-bed interface. This has three important implications for any attempts to compare the model results with reality. First, it is not necessary to form a fracture along the bed, simply to lift it. Second, interactions with the basal drainage system will have potentially large influence on the fate of water arriving via a hydrofracture. Third, enhanced slip at the bed will not necessarily be confined to local 'blisters' as implied by the model.

Each of these three considerations mean that real-world drainage events can play out in ways not simulated in the model. The points about basal drainage and basal slip are especially important, as the authors invoke concern about sea level rise as a justification of their work. In the frozen-bed scenario assumed in the model, any slip will be local and bounded by the surrounding frozen bed; in the real world, slip is less constrained, and can either be enhanced (if the extra water increases areas of ice-bed separation) or reduced (if the extra water encourages development of an efficient conduit). Both of these effects have been described as consequences of surface to bed drainage in Greenland, and show that interactions between hydrofractures and basal drainage are of great importance. 

I am not suggesting that the authors include a basal drainage component in their model. Science often proceeds incrementally, and it is unreasonable to expect that this study should solve the complete problem. I am suggesting though, that the authors more fully acknowledge the complexity of the full problem, and more carefully identify which issues their paper has explored and which ones remain unsolved. It is also worth highlighting that the model also requires the location of initial crack to be specified, a limitation that will need to be overcome before exploring the kind of "flood wave" scenario the authors invoke in lines 495 and forward.

These comments aside, I like this paper, and commend the authors on some interesting work. With the addition of a bit more context (and subject to the technical issues raised by the other referee), I recommend publication.