Review of Trunz et al. 2022

General Comments

Well-written, timely, and based on a substantial body of recent work this manuscript would make an excellent contribution to The Cryosphere and the topic in general. The manuscript presents the results from a moulin-channel model - the methods for which are described in a previous study (Andrews et al., 2022). The model does not include any glacier hydrological systems other than a single moulin and a single Rothlisberger channel and depending upon your perspective this represents the main strength or the main weakness of this study. The conclusion that there is additional “base flow” that contributes to and damps flow within moulin-connected channels is uncontroversial given that many moulins occur on the Greenland Ice Sheet and many appear to be connected (e.g. Andrews et al., 2014). However, this does not take away from the contribution this study makes in very neatly explaining observed moulin head variability using a numerical model. The assertions made regarding differences in “base flow” in the lower and upper ablation areas (Section 5.3.3) are intriguing and demonstrate recent advances in our understanding of the moulin-connected drainage system. 

The main finding of this manuscript is that channel growth is too slow - and channels are therefore too small - to explain the observed damping of diurnal moulin head variations, when the model is forced by local moulin inputs alone. This finding is somewhat similar to that in Dow et al. (2014) but a different interpretation is given in that study: that channels are unlikely to form or persist at high elevations. Current evidence (e.g. Covington et al., 2020; Chandler et al., 2021) as well as the modelling presented in this manuscript suggests channels do form at high elevations and that connectivity to other channels may help them persist. Can the discussion of this be expanded slightly?

Specific Comments

L36 - Water pressure in moulins was also measured by Holmlund & Hooke (1983) and Vieli et al. (2004).

L396 - Section 5.3.2 and Sections 5.3 and 5.4 in general - As mentioned above, a limitation of this study is its application of a moulin/channel model without including other components of the subglacial drainage system. This is stated clearly in Section 5.3 (L375) but isn’t mentioned in Section 5.3.2 which deals only with subglacial channels. It would be clearer if Section 5.3.2 was renamed from “Subglacial network connectivity and base flow” to something more specific to channels e.g. “Subglacial channel network connectivity and base flow”. 

Note that filling of a moulin from the base by subglacial water was observed by Holmund & Hooke (1983) and has been observed in boreholes by several studies (e.g. Gordon et al., 2001) as well as being the focus of Mejia et al. (2021). This suggests that “reverse flow” into moulins does occur and is unlikely to be limited to flow within channels. I appreciate that the model cannot reproduce this - and I don’t suggest this is attempted in this paper - but can this limitation be discussed alongside the direct evidence listed above? 

Section 5.3 appears to overlap with Section 5.4. Is Section 5.3 only concerned with damping caused by storage within the unchannelised system and not damping caused by recharge from the unchannelised system? Can this be made clearer? Can these sections be combined?

In the Section 5.4 heading what does “external” relate to? Is it the same as “non-local”? Does it mean from an unchannelised system.

As in previous studies (Dow et al., 2014; Meierbachtol et al., 2014) Shallow surface and bed slope have a critical role in channel development. How was the slope ratio of 0.01 (equivalent to ~0.6 degrees) measured? Is the bedslope assumed to be the same as the surface slope in Table 1? If so, state this in the methods. If you were to increase the surface and/or bedslope would this change the results?

L128 - Stating the p-value for correlation as a measure of accuracy is inappropriate. The low p-value suggests that the correlation between modelled discharge and measured water level is unlikely to be due to random variation. The p-value does not tell us the degree of accuracy as a statistically significant correlation is plausible for any variables that co-vary regardless of the magnitude (or units) of the variables, or whether there is a causal relationship. For the same reason the statement of “agreement” on L124 between the same variables as above is not strictly speaking supported by the coefficient of determination given, which is a measure of correlation rather than accuracy, though this depends on the intended meaning of “agreement” which is relatively vague. Strictly speaking, the coefficient of determination of 23% suggests 23% of the modelled discharge can be explained by variation in stream level. Reporting the correlation between measured stream level and discharge may be useful but it should not be described as accuracy.

A better measure of accuracy would be the root mean square deviation between modelled and measured discharge. To be clear, I see no problem with the modelled discharge imperfectly matching the measured discharge given how difficult it would be to model rainfall and turbulent heat fluxes. 

L129 - Can you state that m and c are linear regression coefficients and state which variables are the subject of the regression? Is this the melt model calibration mentioned on L123? It’s unclear why day of year 205 was included in the regression even though it was affected by rainfall. Overall, this section needs expanding and revising, possibly with the addition of a figure showing the linear regression. 

L216 - can you revise this sentence? The melt model reproduces melt not “particular weather conditions”. It’s unclear what is being underestimated. If correct, can you revise to make it clearer that melt is underestimated by the model under cloudy conditions and consider adding an appropriate citation.

L230/L251, Fig. 8 and Fig. 8’s caption - consistent description of the error bars would help the reader.

L244 - Is “ranges” required here because you refer to the ranges in amplitude or is range already implicit in amplitude?

Fig. 9 is difficult to interpret. The y-axis label gives two versions of A_h with the modifier in brackets meaning different things (0 relates to amplitude without baseflow, while n relates to time lag). Could A_h(0) be written as A_h minus A_bf? A different symbol for time lag to n, which is usually sample size, would be more intuitive. The symbol for time lag needs to be used consistently (e.g. in the text, legend and caption). A simpler y-axis label could be used and defined in the caption e.g. something like “normalised diurnal head amplitude”. The plot could be labelled with arrows showing where on the x-axis A_bf exceeds A_in and vice versa.

L274 and other occurrences - consider using the past tense for things that were done in the past. 

L367 - This is the orthodox view that crevasses are unlikely to be sufficiently open below a few tens of metres below the surface but there is not that much observational evidence to support this view. It is worth remembering that the moulin would have originated from water flowing into a fracture and that englacial conduits often follow such fractures (e.g. Gully, 2009). I don’t think storage within englacial fractures can be definitively ruled out. See also evidence for fractures at depth presented in Hubbard et al. (2020) and evidence for energy released by refreezing meltwater which potentially occurs in fractures (Luthi et al., 2015).

L451 - reference Figure 10

L452 - discrete in what sense, temporal or spatial or both?

L459 - Nienow et al. (2005) comes to this conclusion by interpreting velocity data on the basis of pressure measurements presented elsewhere (e.g. Hubbard et al., 1995; Gordon et al., 1998). It may be better to cite the studies with hydrological observations of this process.

L484 - Can you reconcile the assertion that at lower altitudes subglacial water flow is steadier while at higher altitudes it is lower in magnitudes but greater in diurnal oscillation amplitude, with the opposing observations in Covington et al. (2020)? Perhaps this would apply if fixed moulin and channel sizes were assumed. Could you expand this assertion in the conclusions slightly to better reflect how it is described (very nicely) on ~L445?

Technical Comments

L17 - Move citation beginning “(Yang …” to before “that”

L74 - delete second citation to Morlighem et al. (2017)

L90 and other occurrences - check whether “Fig. 3” should be “Figure 3” when referenced in the main text outside of brackets. 

Fig. 3 and other occurrences in Figures - units should not be in italics. 

L112 - the first “(Q_p)” is unnecessary. 

L112 - Consider using a different letter for coefficients, or being consistent with the subscripts, or otherwise reducing the potential for confusion; currently there are two C_x’s for concentrations (L95), a C_p for the peak discharge coefficient (L112), and a C for the runoff coefficient (L115). 

L137 - delete “and” and add a comma before phi

L152 - Figure 4 and its description could be easier to follow. In the text, figure, and caption can the same terms be used for “turbulent-underwater” and “open-channel” melt be used? The description of moulin deformation modelled as viscous and elastic needs to be separated from the shear deformation of the ice modelled using Glen’s flow law. Can you mention that the former is modelled using a Maxwell model. Overall the description of the model needs to be expanded to briefly introduce all of the components without the reader needing to refer to Andrews et al. (2022) to make sense of the modelling approach. 

Figure 5 caption - In “(d-e) Moulin shape evolves with surface input.” does this mean that other processes affecting moulin shape were excluded?  

Should “Simulations EMa”be singular? And should “head of oscillation” be “oscillation of head”?

Generally interpretation of figures should be left to the main text and not included in the caption. 

Figure 5 c,e,g y-axis labels - Symbol z has not been defined in the text, and in Figure 4 “Height above the bed (m)” is written out in full.

L193 - specify the sim for the fixed 5 m model run.

L206 - plural ‘radii”

L214 - Specify the sim before referred to.

L254 - Specify Sims E1a-d.

L259 - “of oscillation” could be omitted.

Fig. 8 - Units incorrectly in italics. 

L290 - is it necessary to specify normalized here when discussing in general terms? Is this A^*_in?

Table 3 - Can you use a narrower dash to indicate that a variable is unitless to avoid ambiguity with the minus sign, or perhaps just leave the units for unitless variables blank as on L353? (A narrower unitless symbol is used on Figure 9).

L333 - add “moulin” after “FOXX”

L380 - add “respectively” after “base flow”

L477 - revise to avoid the apparent contradiction in dismissing “subsurface inputs” which would dismiss “basal inputs”.

L468 - the following sentence needs revising to make sense “with only a small portion of the water can” 

Additional References

Gordon S et al. (2001) Borehole drainage and its implications for the investigation of glacier hydrology: experiences from Haut Glacier d’Arolla, Switzerland. Hydrological Processes 15, 797-813. doi:10.1002/hyp.184.

Gulley J (2009) Structural control of englacial conduits in the temperate Matanuska Glacier, Alaska, USA. Journal of Glaciology 55, 681-690.

Holmlund P and Hooke R (1983) High-water pressure events in moulins, Storglaciaren, Sweden. Geografiska Annaler 65A, 19-25.

Hubbard B et al. (2021) Borehole-Based Characterization of Deep Mixed-Mode Crevasses at a Greenlandic Outlet Glacier. AGU Advances 2, e2020AV000291. doi: 10.1029/2020AV000291.

Luthi MP et al. (2015) Heat sources within the Greenland Ice Sheet: dissipation, temperate paleo-firn and cryo-hydrologic warming. The Cryosphere 9, 245-253. doi:10.5194/tc-9-245-2015.

Meierbachtol T, Harper J and Humphrey N (2013) Basal drainage system response to increasing surface melt on the Greenland Ice Sheet. Science 341, 777–779. doi:10.1126/science.1235905.

Vieli A et al. (2004) Short-term velocity variations on Hansbreen, a tidewater glacier in Spitsbergen. Journal of Glaciology 50, 389-398. doi: 10.3189/172756504781829963.

Observed and modeled moulin heads in the Pâkitsoq region of Greenland suggest subglacial channel network effects

Trunz et al.

General comments

This is an interesting and well-written paper that tackles the problem of how to interpret moulin water pressure (head) records. The paper addresses an important topic because moulins are a less inconvenient way of observing subglacial hydrology, than methods such as boreholes.

The main focus of the paper is a suite of experiments with the MouSh model coupled to a subglacial channel model. The experiments attempt to match a moulin head record obtained in the well-studied Paakitsoq region of Greenland. The crux of the problem seems to be creating enough damping in the simulation, so that simulated diurnal head variations remain small enough. A match can be optimised by specifying a suitable subglacial base flow or a very large moulin shaft. Although the experiments themselves are interesting conceptually, there are two very significant weaknesses that limit how useful this study is in its present form.

First, to match the simulated moulin head with observations, the authors need to add subglacial water (base flow). Although this is of course quite reasonable when modelling a moulin that joins a wider scale drainage system, it’s not clear how much of the discrepancy between modelled/observed moulin head is due to (i) the subglacial base flow and subglacial channel model, (ii) uncertainties in the moulin model, and (iii) the rather poorly prescribed melt input. Taking these in turn:

In the case of base flow, it is confusing that it is prescribed as a flow rate (m3/s) even though that is on first impression an obvious choice. The moulin head is a proxy for subglacial water pressure (Pw), not flow, and even though these quantities are related, prescribing flow rather than Pw presumably requires some speculation of drainage system characteristics to calculate Pw in an extra step. In fact without any base flow or stream input at all, the moulin head could simply reflect changes in Pw driven entirely non-locally, provided the moulin remains hydrologically connected. This aspect needs to be clearer so we know what assumptions are needed and how are the associated parameters constrained. It is also not obvious how Qout is calculated at the bottom of the moulin, nor how the subglacial channel model is driven in the context of wider-scale drainage evolution through the season.

A minor related point: the upper limit for baseflow of 5 m3/s in the experiments seems very low for a typical Greenland catchment or even a small part of one.

For the moulin model we need to remember that MouSh is effectively unvalidated, even though it is quite a detailed model. Inevitably there are many unknowns here – for example, my understanding is that MouSh assumes an initial moulin that is a vertical cylinder, which I suspect is far from reality given that most moulins appear to form as hydrofractures in Greenland, which are initially planar. Anybody that has tried lowering sensors into a moulin can guess they are not vertical shafts! Limited exploration confirms that. Perhaps several moulins are connected englacially within an initially planar fracture zone, before even reaching a subglacial channel. Almost like a ‘distributed’ englacial drainage system that might also produce the dampened head record sought by the authors. Of course this is speculative, maybe moulins are vertical shafts after all, but in the experiment we should account for our lack of knowledge in this respect.

Melt input: this is the last important source of uncertainty, and it seems very much brushed under the carpet in Section 3.1.1. Even the very short (2-day) observed time series is not well represented by either the ‘Modeled’ or ‘Idealized’ time series (Fig 3), and the extended parts of the Q time series shown in Fig 3 do not follow the trend of the steam water level, which is discouraging. The vast majority of the melt input time series is extrapolated outside of this short (poorly matched) tuning window. I strongly disagree with their “good confidence in our derived runoff values R and thus our model forcing Qin,model”. I also wonder why most simulations (Tab 2) use the modelled Qin, even though it looks like a worse fit than the idealized? This poor match is not necessarily a problem when constraining a model with observations, if it is clearly acknowledged and if the uncertainty it introduces is more rigorously accounted for.

I think to address this first weakness we would need some detailed error analysis before any comparison between the simulations and observed head time series can be interpreted in a useful manner. I would envisage an ensemble (e.g. Monte Carlo or latin hypercube sampling) or more advanced statistical approach to account for the very many uncertainties in the MouSh model, the treatment of baseflow, and the stream inputs. The long discussion section (which is currently very speculative without even a basic error analysis) should then focus on how much of the record can be confidently interpreted in terms of real variations in subglacial water pressure (or drainage characteristics), and how much cannot be untangled from uncertainty in the simulation and inputs.

The second weakness is its relevance, which of course I acknowledge is limited by where field data are available. In Fig 5 it is evident that the study is conducted mid melt season when there is relatively little variation in moulin head (range looks like 250 to 400m, but much of period close to 300m), and as such Pw is always well below the ice overburden pressure. In fact the authors choose a period mid-season when h is around 60% of H and varies diurnally by about 10%. In these conditions we would expect the moulin head to have a minimal effect on ice motion. Probably there are some data for this region that could answer that more precisely. What the community needs more, I would suggest, is for the study to simulate the early season formation of moulins as part of efforts to simulate the duration/extent of ice acceleration in spring. What role do moulins play in the evolution to channelised drainage & associated ice deceleration? It’s not clear how the results would help that aim, at the moment. However, some simulations of moulin head very early in the season could provide valuable pointers for interpreting moulin pressure records in the more dynamically important part of the cycle, or in thicker ice. Similarly, experiments simulating extreme melt pulses mid or late season could be useful. But, related to the first point, we would need to know how the subglacial channel model is driven by / coupled with the larger-scale hydrological evolution.

Overall I think this could go one of two ways – (i) keep the focus on the link with available field observations, by carrying out some detailed error analysis, or (ii) accept that the field observed melt inputs are perhaps too limited/uncertain at present, and instead focus on a conceptual study that is not tied to that one location and can explore controls on head variations across a wider (more interesting) range of sites/conditions. Either of these directions could turn this into a really useful study to help interpret or design experiments with moulin water pressure records.

Minor comments from the introduction

L2: I believe water pressure also influences motion in some marine terminating glaciers (not just land terminating).

L19 I think some observation based papers showed the influence of temporal melt inputs, before Schoof 2010 – worth to cite these here too.

L28 again there are earlier papers describing drainage evolution in Greenland (Bartholemew et al 2011 EPSL?), and of course even earlier elsewhere (1980s/90s work on alpine glaciers).

In general there is a lot of self citation in the intro that could be expanded to include earlier work from other groups.

Fig 8 confusing that red dots can be either observations or simulations. Can sims not be blue or some other colour?