the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Ice viscosity governs hydraulic fracture causing rapid drainage of supraglacial lakes
Abstract. Full thickness crevasses can transport water from the glacier surface to the bedrock where high water pressures can open kilometrelong cracks along the basal interface, which can accelerate glacier flow. We present a first computational modelling study that describes timedependent fracture propagation in an idealised glacier causing rapid supraglacial lake drainage. A novel twoscale numerical method is developed to capture the elastic and viscoplastic deformations of ice along with crevasse propagation. The fluidconserving thermohydromechanical model incorporates turbulent fluid flow and accounts for melting/refreezing in fractures. Applying this model to observational data from a 2008 rapid lake drainage event indicates that viscous deformation exerts a much stronger control on hydrofracture propagation compared to thermal effects. This finding contradicts the conventional assumption that elastic deformation is adequate to describe fracture propagation in glaciers over short timescales (minutes to several hours) and instead demonstrates that viscous deformation must be considered to reproduce observations of lake drainage rate and local ice surface elevation change. As supraglacial lakes continue expanding inland and as Greenland Ice Sheet temperatures become warmer than 8 °C, our results suggest rapid lake drainages are likely to occur without refreezing, which has implications for the rate of sea level rise.
 Preprint
(3076 KB)  Metadata XML

Supplement
(3690 KB)  BibTeX
 EndNote
Status: closed

RC1: 'Comment on egusphere2024346', Anonymous Referee #1, 16 Apr 2024
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 ‘endmember’ 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 nontrivial  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 nonsmooth, 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 sticksliplike 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 waterlevel 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 preexisting 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 zerofracturetoughness 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 icebedrock 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 latertime 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 nonlinear viscoelastic 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 nonNewtonian 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 nonlinear viscous deformation (nonrecoverable, and given by a shearthinning model). A ‘viscoplastic’ model would typically be taken to indicate that there is a plastic ‘yield’ stress, below which there is no nonrecoverable deformation; at which the material deforms plastically; and above which the material flows viscously (see e.g. much literature on viscoelastoplastic 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 timescales’ which could mean anything.
 p.6 equation (3) has some weird typesetting on the third line (missing equals?)
 p.9 The authors choose a ManningStrickler 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 laminartoturbulent 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 linearelastic 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 (1618) 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 (tt_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?
Citation: https://doi.org/10.5194/egusphere2024346RC1 
AC1: 'Reply on RC1', Emilio MartinezPaneda, 25 May 2024
The comment was uploaded in the form of a supplement: https://egusphere.copernicus.org/preprints/2024/egusphere2024346/egusphere2024346AC1supplement.pdf

AC1: 'Reply on RC1', Emilio MartinezPaneda, 25 May 2024

RC2: 'Comment on egusphere2024346', Douglas Benn, 02 May 2024
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 'viscoplastic' is not appropriate and should be changed to "viscous" or "nonlinear viscous" throughout. It is also unnecessary to refer to the "socalled 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 "socalled" 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 icebed 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 realworld 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 frozenbed 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 icebed 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.
Citation: https://doi.org/10.5194/egusphere2024346RC2 
AC2: 'Reply on RC2', Emilio MartinezPaneda, 25 May 2024
The comment was uploaded in the form of a supplement: https://egusphere.copernicus.org/preprints/2024/egusphere2024346/egusphere2024346AC2supplement.pdf

AC2: 'Reply on RC2', Emilio MartinezPaneda, 25 May 2024
Status: closed

RC1: 'Comment on egusphere2024346', Anonymous Referee #1, 16 Apr 2024
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 ‘endmember’ 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 nontrivial  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 nonsmooth, 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 sticksliplike 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 waterlevel 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 preexisting 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 zerofracturetoughness 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 icebedrock 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 latertime 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 nonlinear viscoelastic 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 nonNewtonian 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 nonlinear viscous deformation (nonrecoverable, and given by a shearthinning model). A ‘viscoplastic’ model would typically be taken to indicate that there is a plastic ‘yield’ stress, below which there is no nonrecoverable deformation; at which the material deforms plastically; and above which the material flows viscously (see e.g. much literature on viscoelastoplastic 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 timescales’ which could mean anything.
 p.6 equation (3) has some weird typesetting on the third line (missing equals?)
 p.9 The authors choose a ManningStrickler 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 laminartoturbulent 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 linearelastic 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 (1618) 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 (tt_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?
Citation: https://doi.org/10.5194/egusphere2024346RC1 
AC1: 'Reply on RC1', Emilio MartinezPaneda, 25 May 2024
The comment was uploaded in the form of a supplement: https://egusphere.copernicus.org/preprints/2024/egusphere2024346/egusphere2024346AC1supplement.pdf

AC1: 'Reply on RC1', Emilio MartinezPaneda, 25 May 2024

RC2: 'Comment on egusphere2024346', Douglas Benn, 02 May 2024
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 'viscoplastic' is not appropriate and should be changed to "viscous" or "nonlinear viscous" throughout. It is also unnecessary to refer to the "socalled 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 "socalled" 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 icebed 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 realworld 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 frozenbed 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 icebed 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.
Citation: https://doi.org/10.5194/egusphere2024346RC2 
AC2: 'Reply on RC2', Emilio MartinezPaneda, 25 May 2024
The comment was uploaded in the form of a supplement: https://egusphere.copernicus.org/preprints/2024/egusphere2024346/egusphere2024346AC2supplement.pdf

AC2: 'Reply on RC2', Emilio MartinezPaneda, 25 May 2024
Viewed
HTML  XML  Total  Supplement  BibTeX  EndNote  

363  110  26  499  48  31  16 
 HTML: 363
 PDF: 110
 XML: 26
 Total: 499
 Supplement: 48
 BibTeX: 31
 EndNote: 16
Viewed (geographical distribution)
Country  #  Views  % 

Total:  0 
HTML:  0 
PDF:  0 
XML:  0 
 1