This paper introduces a new model, MouSh, designed to emulate the seasonal evolution of moulins and examine the impact of this on the basal drainage system. A model like this has been long overdue and it’s exciting to see it implemented to investigate both the evolution of the englacial system but also the subglacial system, which may be controlled strongly by moulin capacity and fill and drain rates.

The paper is nicely written and details many aspects of the model which has a lot of components. Looking at the results for the viscous deformation, turbulent melting and the application of Glen’s Flow law provides interesting information about the evolution of these features and how they might significantly change radius even on a daily basis, which is not something that I was aware could happen.

I do have some major concerns which I list below but, once these are taken into account, I think this will be a strong addition to glaciological modeling.

Major issues

1. I’m confused by the implementation of elastic equations in this context. Most applications of elastic equations in glaciology, that I’m familiar with, are for a situation with a bending beam or plate (e.g. ice shelf flexure) in response to a changing force. I have a difficult time understanding how this can apply from constant ice force inwards into a moulin in an elastic rather than viscous form, particularly since elastic deformation should be instantaneous with changing force but here it is from the change in resistance (the water), and I’m not convinced that they are equivalent. As far as I know, this type of calculation with both elastic and viscous deformation in a moulin or borehole is new to glaciology and so the approach needs more explanation/justification.

Part of the problem is that equation 4 is difficult to follow in its current form. If this is a new way to apply elasticity to a moulin then the equation needs to be fully derived in an appendix. If it has been applied in glaciology before you need citations. You base your equations from Aadnoy (1987), but this isn’t even included in the reference list.

Both the elastic and viscous deformation rates that you plot are higher than I would have assumed for this situation. Previous analysis on borehole deformation has closure rates one or two orders of magnitude smaller (e.g. Paterson, 1977). Also Catania and Newmann (2010) argued closure would primarily occur in the base of the moulin and not the top 80%. A discussion on why yours are so much higher and/or examples of other systems that deform this rapidly would help.

2. Following on from this, where have you got your Maxwell times of 10-100 hours from? I believe the Maxwell time should be more in the range of a few hours. Therefore viscous deformation should be the primary application for moulin shape evolution. I’m also very unclear how you’re transitioning between elastic and viscous deformation with this model and applying the Maxwell time. It seems these are both being calculated separately but continuously given that you are plotting both over sub-hourly, multi-day timescales?

3. The treatment of turbulent melting and refreezing is confusing. Why only include refreezing outside of the melt season? Do you assume no refreezing overnight during the melt season? I see you say refreezing occurs only when water flow is laminar but it’s not clear to me that water flow will be always turbulent from the beginning to the end of the melt season. This needs more justification in the text by reporting the expected Reynolds numbers.

4. I understand why you’ve applied a simplified subglacial model given the complexity of the moulin model. However, both the description of the basal channel model and the application are confusing. From what I can gather you’re calculating channel characteristics at the moulin outlet (using ice pressure and moulin head pressure) but are applying a constant hydraulic potential gradient from the moulin to the terminus, so only producing one output point. The length scale calculations from moulin to terminus are not ideal in application to a continuously evolving channel (which will not be linear in terms of pressure) and are likely unrealistic. Instead why not apply a range of hydraulic potential gradients to test how those impact the moulin evolution? That would be much clearer to show how the pressure change at the bed impacts and is impacted by the moulin head.

One of the significant concerns I have about the channel is the necessity of a large base flow. Looking at Figure 7 the base flow is the main driver for channel evolution and at an input rate of ~20m3/s that’s not surprising it’s the primary control. To better determine the role of the changing moulin head it would be better to avoid adding additional time-varying water inputs at the bed since it’s not clear that it’s at all realistic. Instead, a static background water flux and/or a larger initial channel size could help with stability issues.

It’s generally hard to believe the channel outputs that you present as it’s not clear what the differences are between basins and experiments in terms of the hydraulic potential gradient, and because of that large base flow rate. However, as this paper focuses on the moulin model, so should the results and discussion. The role of a basal channel in this case is to present semi-realistic evolution characteristics to feedback with the moulin water levels. This does not give you much information about what is happening at the bed anywhere downstream of the moulin so that should not be widely discussed. Along these lines, before you begin your moulin model methods, you look at subglacial channel routing in section 2.1, which is misleading for the reader. This section does not seem relevant to this paper because of the highly simplified nature of the channel model that you apply and it would be better to start the paper with the moulin model methods.

5. The discussion at the moment focuses a lot on how moulins are formed, the subglacial system, and englacial void ratios. These don’t seem directly relevant to your main findings from this complex model, which are the changes in shape, melt rate and deformation of the moulin. Particularly given that the subglacial model is much more simple and this is the first step in coupling to a more dynamic subglacial model, the discussion in this paper should be focused on the moulin evolution. There are many interesting outputs from your model runs that you could discuss in terms of the deformation of the  ice possibly moving the input of the moulin at the bed along with stretching the length of the moulin; where in the moulin and at what time of the season water would be stored at higher or lower pressure influencing the subglacial system; the influence of the moulin shape on the head etc.

Line-by-line comments

11. where does that 10-14% number come from?
12. they constitute most of the englacial system – what about englacial channels?

      38-39. you already said this in your first sentence of the introduction.

71. what do you mean ‘relative path length’? This whole paragraph is confusing because you’re discussing basal hydraulic potential rather than moulins.

      84-95. I’m very confused. what do you mean by theoretical flow accumulation? Are you saying you’re defining the catchments at the base of the ice? You’re defining channel lengths at the bed? What is a subglacial channel node? Channels should join up dendritically towards the terminus in any case and are therefore linked rather than in separate segments. 

104. why do you initiate with a semi-circular, semi-elliptical shape? There doesn’t seem to be any reason for this and it primarily serves to complicate your equations and your analysis.
105. undefined parameters for Maxwell time. Where do you get the equation in parenthesis from?
106. if moulins form by drainage into crevasses and hydrofracture why do you assume it’s compressive?
107. you say viscous deformation is the dominant process over a 1 day timescale but you plot your viscous deformation on much smaller timescales showing diurnal variation.

           154-156. specify here this opening and closure is relative to the pressure difference at depth – moulin should not open at all depths when above flotation – only in regions of the moulin where the relative pressure is higher than flotation. Looks like you’re calculating this in the next section but this should be clarified here.

160. why have both equation 5 and equation 6?
161. laminar flow is when the Reynold’s number is less than 2300.
162. do you mean all ice is at the melting point, not water?
163. what difference do you find with these alternative approaches?
164. you haven’t told us about the equation you use for turbulent melt. I see that you have it later in the paragraph but it’s confusing in this line because it implies we already know how you calculate it.
165. but you do assume it’s at the pressure melting temperature in your refreezing section. I’m getting confused.
166. how modest?
167. what is S in equation 18? Again I thought the ice on the moulin wall was at the pressure melting point?
168. presumably the unit hydrograph is to allow a lag for the runoff to reach the moulin? If so, you should state that.

    322-329. this last paragraph seems more appropriate in the next section

325. specify what you mean by englacial void ratio here? Why would that impact flow from upstream?
326. what two elements?
327. assuming b is elevation above sea level, why include it if you have zero bed slope?
328. where is the base flow added?
329. you mentioned Qin above in section 2.2.4.1, which included the baseflow. That’s not being added directly into the moulin I assume? This should be clarified.
330. You need a justification for your choice of enhancement factor. It seems like you’re applying this factor between 1 and 9, but measurements by Luthi et al (2002) in Greenland ice suggested it can reach up to 2.5 in Holocene ice but is closer to 1 above that depth.
331. what are these basins?
332. why are these lengths so different? I don’t think these lengths scales help your model application – it makes it confusing (see comment at the beginning).
333. how can they reach equilibrium with constantly changing water input and a constantly evolving basal channel? Or is this with a constant input?

     440-445. Changing Youngs modulus for elastic expansion increased moulin volume by 38% and capacity by 56%? This seems like a much more significant change than I would assume from elastic deformation in this context.

461. it would be more appropriate here to say the outflow is sensitive to the steepness of the basal pressure gradient.
462. but the increase in water level should increase the pressure gradient and cause faster flow through the subglacial channel and melt opening?
463. this discussion and Figure 6 show a diurnal change in viscous deformation by up to 20cm and 10cm elastically. Then the diurnal phase change up to 30cm/day. Is that saying you melt up to an extra 30cm a day? And that the moulin diameter pulses in and out every day due to melt countered by viscous and elastic deformation? Is there any evidence for this from moulin measurements.
464. what would cause a moulin to change radius diurnally more in thicker vs. thinner ice?
465. the similarity between Basin 1 and 2 is likely because the basal channel model is driven by a similarly large background influx rather than by the changing conditions of the moulin itself.
466. how could thick ice viscously close channels if water is above overburden pressure?
467. specify which system sees the increase.
468. you have elastic processes in the channel too? Any references to show this is justified in basal channels?
469. I’m unsure why you’re discussing initial moulin formation processes which aren’t the focus of your study. Moulin evolution, yes, and you have plenty of interesting things to talk about on that subject.
470. exploration would be good to validate your model so I wouldn’t discount it.
471. can you clarify this sentence. You are saying field measurements show 103-112% of overburden, or 3-12% of overburden? The former seems more likely. But 20% above overburden? That should be on the surface?
472. some boreholes hit more efficient systems, as explored by Meirbachtol et al (2013).
473. rephrase ‘variations in diurnal water level variability’
474. I don’t think an englacial void ratio is used to resolve diurnal basal pressure. How would you get a spatially variable englacial void ratio? What has this got to do with moulins?
475. you said earlier in the paper that moulins are used as source inputs for models. How does this link to englacial void ratio? The change in water level is a moulin because of increase/decrease in diameter will impact the water supply to the base via the pressure. Perhaps you mean a storage parameter in models? I certainly think it’s worth coupling with subglacial hydrology, but I’m not sure this paragraph makes sense. Line 696 covers this possibility and is an important point to make.
476. what do you mean a static shape instead of static cylinder? The Trunz et al, in review paper is mentioned a lot which is frustrating since we don’t have access to see what it discusses.

     731-739. see my above comments about elastic deformation. You need more justification for these statements given that it hasn’t been included in subglacial models to date.

742. are you sure it’s not that subglacial channels form where there are moulin inputs?
743. rephrase this sentence.

Figures

Figures 4 and 5. Why does the y axis of the diurnal range go up to 0.4 if values don’t go above 0.2?

Figure 6.  I’m intrigued by the shapes in f. Why is there more turbulent melting in the middle of the borehole? What are the factors contributing to the differences between elastic and viscous deformation shapes and rates? In g since it seems to have reached equilibrium within a day or so it would be useful to zoom in so we can see the lines better.

Figure 7. In your thickest ice example for moulin water level, it looks like your moulin is overflowing. Also in d) where is the 741m example? The channel size looks almost entirely dictated by the background flow you input with small diurnal variability on top.

Figure 8. This is a really interesting figure. Why not discuss the shape changes (particularly due to Glen’s Flow law) more in the manuscript?

Figure 9 b. What happened around day 32?

References

Catania, G. A., & Neumann, T. A. (2010). Persistent englacial drainage features in the Greenland Ice Sheet. Geophysical Research Letters, 37(2).

Lüthi, M., Funk, M., Iken, A., Gogineni, S., & Truffer, M. (2002). Mechanisms of fast flow in Jakobshavn Isbræ, West Greenland: Part III. Measurements of ice deformation, temperature and cross-borehole conductivity in boreholes to the bedrock. Journal of Glaciology, 48(162), 369-385.

Meierbachtol, T., Harper, J., & Humphrey, N. (2013). Basal drainage system response to increasing surface melt on the Greenland ice sheet. Science, 341(6147), 777-779.

Paterson, W. S. B. (1977). Secondary and tertiary creep of glacier ice as measured by borehole closure rates. Reviews of Geophysics, 15(1), 47-55.

Apologies for taking a long time to produce this review.

This is a very interesting paper that presents a new model for the evolution of moulin geometry, and explores how the results of this model for moulin shape and water level depend on various model parameters.   It is argued that moulins comprise a sizeable fraction of the englacial-subglacial drainage system in Greenland, and that the time-evolution of their volume is a potentially important feature to include in englacial/subglacial models, offering improvements over a model that assumes a static moulin volume.

The study is an interesting one and I believe it deserves publishing in some form.  However, I do have quite a lot of detailed questions, and some concerns, about the ingredients that go into the model.  I will focus this review largely on these model details, from section 2 of the paper.   Some of these may be sorted out by clarification as to what equations have actually been solved.  As a general comment, there appears to be quite a lot of duplication of notation, which overcomplicates the presentation of the model and causes some confusion.  I think it would also be helpful to express the physics in terms of differential equations rather than discrete increments that implicitly include time-steps.

Major comments

Section 2.2.1 - the rationale for modelling the moulin cross-section with this strange egg shape was weak for me.  It significantly complicates the model to do this, rather than to assume it has a circular cross-section, and it was not at all clear to me that there was any great advantage in doing so.  It is also not clear how r_1 and r_2 are separately evolved, and this needs to be made clearer.  I ended up with the impression that the difference is likely because the open channel flow above the water line gives rise to a change in one of these but not the other; but below the water-line it seemed that r_1 and r_2 would evolve identically and therefore stay the same, assuming they start the same?  However, this should be made clear by telling us what exactly are the equations that govern the evolution of r_1 and r_2.  I would, at the same time, encourage the authors to think about simplifying things and assuming circular symmetry, since I think many of the results would still apply, and I think it would give a model that is more likely to be adopted by others.

Section 2.2.2 - I had great difficulty following the treatment of elastic deformation, and am slightly concerned that this is not dealt with correctly.  In particular, a number of figures (figure 6, figure 9) compare viscous and elastic ‘deformation’ as a *rate*, with units m/d.  Elastic deformation is not a rate - it is an instantaneous deformation and it results in a displacement (relative to some reference state) that is fixed, for fixed stress - in this context, that is the change in radius given by (4).  Presumably this must be viewed as relative to some ‘reference’ radius that can evolve in time due to viscous deformation and phase change.  The elastic displacement in (4) does itself evolve in time due to changes in water pressure, and therefore gives rise to a deformation rate that is d(\Delta r_E)/dP * dP/dt, i.e. proportional to the rate of change of water level, and perhaps that is what is being plotted in these figures, but I did not really have this impression. If that is indeed what is meant, note that the elastic deformation rate depends on the rate of change of P, not on P itself, so whether the pressure is above or below overburden is irrelevant to the sign of the deformation rate (it is instead a question of whether P is increasing or decreasing).

This concern is tied up with the question above of how exactly r_1 and r_2 are evolved.  It seems to me that you would want to have ‘reference’ values of these that evolve according to the viscous processes; they satisfy an equation of the form dr/dt = melt-back - viscous closure (very similar to the subglacial channel in (28)); and then you want to add the elastic deformation given by (4) on top of those evolving reference values to get the actual radius at any instant in time.

In equation (4), I would be inclined to simply ignore the deviatoric stresses, which I expect are relatively small in most cases compared to the effective pressure P (it was not clear to me what you have actually assumed for them in the examples).  Given that you are comparing with a null model which contains no moulin physics whatsoever, I think there is some advantage in not making this one overly complicated!   Note that there are in any case some missing brackets in this equation.  In equation (5), the first term is presumably set to zero for z larger than h_w (i.e. above the water line)?  This equation could be made more consistent with (6), which is essentially the same thing, but where P is now called sigma_z.  (6b) should again by zero for z larger than h_w, I think.

Equation (9) is a strange way of discretising the time-derivative and this is where confusion starts to arise as to how r is actually evolved, because this gives an incremental change in r (both r_1 and r_2 ?) due to only viscous processes, and it is not clear how this is combined with the changes due to phase change and elastic deformation.  The viscous closure of a moulin due to (7) is essentially identical to that for a subglacial channel as described in (27) and as described by Nye (1953) for the closure of a borehole.  I think it would be helpful to express it as a contribution to the time-derivative dr/dt, as (effectively) done in (27).

Section 2.2.2.2.2 (I don’t think I’ve ever seen quite so many subsections!) - The downstream deformation of the ice is interesting, but it wasn’t clear to me how it is incorporated into the model.  It seems like it translates the ‘centreline’ of the moulin?  But doesn’t affect r_1 and r_2? So does it actually have any effect on the rest of the model or is it just relevant for the visualisations like in figure 8?  The formula in (10) assumes no slip at the bed, which is presumably not always going to be the case?

Section 2.2.3 - I was a bit confused why melting and refreezing are treated separately - you could simply write down an energy balance that allows for either to happen automatically, depending on the relative magnitude of turbulent heating and the conduction into the ice, without having to have any ‘switch’ between melt season and not.  In (11), I would have thought that the dT/dx should really be a dT/dr, i.e. the radial temperature gradient away from the (roughly) cylindrical moulin; the distinction between them is quite important because conduction around a point source in two dimensions (ie. in the x,y plane) is very different from conduction in one dimension (i.e. in x alone).  That said, solving the heat equation in the ice for each different z seems a lot of work for a model of a single moulin, and I wonder if a reasonable approach would be to simply *estimate* the temperature gradient at the moulin wall, dT/dr, as (\Delta T)/r_m, where \Delta T is the temperature difference to the far-field ice and r_m is the moulin radius.  That would be consistent with the way you incorporate the estimate of sensible heat in (18) when considering melting.

Is \Delta r_t in (19) the same as the melt rate m in (14)? And what exactly is Q here, in relation to the other Qs mentioned later (Q_in, Q_out, Q_base)?  If I understand the picture correctly I think it ought to be Q_out - Q_base, since that’s the flow out of the moulin into the subglacial channel.  This could all be made clearer with more consistent notation.  I couldn’t follow what is used for the melting in the open channel zone on L287-295; it says you use (17), but that doesn’t seem helpful. I would have thought you want to use something more like (19), but with Q replaced by Q_in, and with a modified hydraulic radius and perimeter.

Section 2.2.4 - Equation (22) needs to include Q_base, similarly to equation (24).  In fact, there seems to be some inconsistency and duplication between (22), (24) and (29).  These equations are all expressing mass conservation, and (22) and (29) are really the same equation (I assume that the m in (29) must include the freezing rate -delta as well).  But I think they should include Q_base if you’re going to include Q_base in (24).  And I think (24) should really include some terms to account for the rate of change of the cross-sectional area (it comes from inserting V_m as the integral of A_m from 0 to h in (29)).

Section 2.2.4.2 - I think it would help to have a schematic picture of the moulin and the subglacial channel showing some of the various variables.  It is slightly frustrating - but I can see that it may be unavoidable - to have the moulin shape model coupled so tightly to a subglacial channel model; ideally you’d like to be able to model the moulin separately.  In this case, it seems that the subglacial channel is assumed to run from the bottom of the moulin to the ice-sheet margin, along which length the channel cross-section would presumably vary in reality, but I think that you assume a single value of S (the value at the bottom of the moulin?) is sufficient to describe how the flow evolves?  This seems a reasonable simplification here, but I think could be explained a bit better, and as I say, a diagram might help.  The ‘b’ in the hydraulic gradient on line 339 seems to disappear when this term is inserted in (28).   The diagram might also help to explain Q_base, Q_in and Q_out.  The use of Q_base seems fine to me, as for most moulins there will likely be water arriving at the bottom of the moulin from upstream as well as via the moulin.

Section 3 - The results section focuses a lot on parameter sensitivity, and it is great that this has been explored so thoroughly, but I found this hard to follow without it having first been outlined some of the general behaviour of the model.  In particular, I think it would be helpful to see some sort of figure showing the periodic states to which the moulin apparently evolves.  Just the fact that the modelled moulin approaches an ‘equilibrium’ does not seem an obvious result, and I think that equilibrium could be described a bit more fully.  Presumably it involves the water level moving up and down on a diurnal timescale, and the moulin opening and closing?   It would also be useful to know how this depends on the moulin input Q_in (for me that would seem more of interest than dependence on drag parameters etc, which we don’t know very well).  It seems quite surprising to me that if such an equilibrium is really reached, it depends on the initial moulin radius.  Also, has the moulin model been run over the course of multiple years (with melt season and a winter), and how does it behave?  This has implications for what an appropriate ‘initial’ moulin radius is, presumably.   I think it would be helpful to have some general discussion along these lines, and figure(s) (perhaps like figure 6 or 8) that show the general behaviour of the model, before going into detail about how certain outputs depend on the parameters, since it would help give those more context.

Figure 6 - see my earlier comments about comparing elastic and viscous deformation.  I just don’t understand what is actually plotted in panel f and g.  Could you express whatever quantity is being plotted in terms of variables in the equations?  Similarly for figure 9, and the associated discussion in section 4.5

More minor comments

L82 - why does taking k = 1 approximate likely channelized pathways?  The usual thinking is that channels would tend to *lower* the water pressure and would therefore be associated with a lower value of k, if anything.

Figure 1 is very nice.  It might be noted that the elastic deformation here is quite different from all the other ones, in that the others are all *rates* - they accumulate every timestep to give continued deformation - whereas the elastic one is just static. 

L185 - the small component of melting due to temperature differences between the water and ice seems to be ignored in the model, since it is later assumed that the water is at the melting temperature ?

L255 - you seem to use both hydraulic diameter D_h and hydraulic radius R_h and it would keep the notation simpler to just work with one or the other.

L262 - it sounds like in the end you take f_R to be fixed (and vary it’s value) so I wasn’t sure what the point of introducing (16) was.

In (18) presumably S is really A_m, the moulin cross-sectional area?

In (19) is dh_L/dz the same as dh_L/dL in (15), and is there significance in the change from lower case to upper case subscripts? 

In (23), time appears to be in hours, not days.

In (24), h is the same as h_w ?

Figure 7 - should there be a purple line in panel (d)?

L662 - I wasn’t able to see this statement about the fixed moulin frequently overtopping the moulin in Fig 11a.  How does the figure show this?