We thank Referee #1 for their comments on our paper titled ‘A new model for supraglacial hydrology evolution and drainage for the Greenland ice sheet (SHED v1.0)’. The comments were constructive and have contributed in making this paper better. In our response below, ‘RC 1’ stands for Comment by Referee #1 and ‘AC’ stands for Author Comments.

RC 1: There are several ways to perform flow routing on a DEM (e.g., D8, Din, multislope etc). Which are you using and how do results vary if you change these? How do you account for DEM uncertainty?

AC: In our analyses, meltwater runoff was routed between the DEM cells, based on the principle of lowest lying neighbour according to which runoff from one DEM cell (i.e., origin) is transferred to a DEM cell that a) lies in the 3-by-3 neighbourhood of the origin and b) has the lowest surface elevation in the neighbourhood.

We also tested flow routing using the standard D8 algorithm according to which, runoff was transferred from the origin to a DEM cell that a) lies in the 3-by-3 neighbourhood of the origin and b) has the highest surface slope with respect to the origin. The results were not substantially different from our approach, and an additional advantage of our approach is that it is slightly more computationally efficient.

We have added text to explain about how we estimate flowpaths between a DEM cell and its Potential Destination Cell (PDC) in section 1 of Appendices and referenced the appendices in the main body of the revised manuscript.

The average uncertainty in elevation associated with the ArcticDEM product that we use in this study is ~0.1 m (Candela et al., 2017). Due to a paucity of in-situ elevation data for our region, we assume 0.1 m to be the vertical accuracy for the DEM of our study area.

We have added a sentence to the revised manuscript giving this uncertainty and have added the citation  Candela et al. (2017) to the reference list.

RC 1: Do you predict evolution of the underlying topography itself? Catchment boundaries, and supraglacial stream networks, likely evolve through the melt season. At the wavelengths of catchments, ice topography may be largely dictated by underlying bedrock topography (e.g., Gudmundsson , connected to stream networks explicitly by Crozier et al 2018). Thermal erosion in streams is rapid and thus dynamic reorganization of the supraglacial network (e.g., stream capture events and other significant/sudden changes to upstream drainage area at a point) are possible (Karlstrom and Yang, 2016).

AC: We simulate the evolution of the underlying surface topography in our modelled channels,by updating the surface elevation of the DEM cells (at every time step) as they are eroded in the ‘Lake overtopping and drainage module’ (i.e., Module 3). Thus it is possible for our model to simulate reorganization of the supraglacial meltwater channels when lakes drain laterally (e.g., Kalstrom and Yang, 2016). However, we do not account for change in surface elevation due to surface melting and ice-flow because we are not currently coupled with an ice flow model, but this would naturally be something we will investigate once we have a coupled model (which is a future goal beyond the scope of the current study). We have now added two sections of text to the revised manuscript to clarify this. This includes an explantion about how we simulate the re-organisation of meltwater channels as a result of rapid lateral lake drainage, as well as an explanation about our plan to move towards coupling with an ice-flow model in the future, so that we can simulate other observed processes (such as re-organisation of surface meltwater channels as a result of rapid vertical lake drainage events). We have also added Karlstrom and Yang (2016) to the reference list.

RC 1: Unchannelized flow through snow and firn is driven by hydraulic head above an impermeable substrate is an unconfined aquifer. If channels are known, one can calculate transit times based on the Dupuit approximation (e.g., Yang et al., 2018). It seems that porous flow through porous ice (not just snow) could be better accounted for – this represents a much slower transport than sheet flow implied by Manning. I was a bit surprised that you did not compare with the model of Yang et al., actually.

AC: When MAR predicts snow cover in a grid cell in our model, flow velocity and travel time is simulated as per Colbeck (1973). The travel time consists of two components. In the first, we estimate the time meltwater runoff takes to trickle down to the base of the snowpack, and in the second we estimate the time taken by the runoff to flow between DEM cells, at the snowpack’s base i.e., assuming saturated flow over impermeable substrate (by using Manning’s equation) (e.g., Arnold et al., 2020). We appreciate that weathered ice may also be porous, and that we do not account for this. However, given that the properties and extent of weathered ice are not known, we think that this is a reasonable approach, which has also been followed by previous studies (e.g., Koziol et al., 2017).

We note that our estimates of modelled flow velocity of meltwater runoff in supraglacial meltwater channels is in the range of 0.001 ms^-1- 0.462 ms^-1, in keeping with values suggested by Yang et al., 2018 (0.2-9.4 ms^-1). The reason our estimates are at the lower end of this range is likely because our estimates are modelled at a spatial resolution of 100 m, and those of Smith et al. (2015), which are used in Yang et al. (2018) are point measurements. We have added text to the revised manuscript stating our average flow speeds in comparison to Yang et al. (2018). We have also added Yang et al. (2018) to the reference list.

RC 1: Does variation of stress matter for assessing hydrofracture initation/propagation? This approach would seem to ignore contributions from bedrock/surface topography? Many depressions where lakes form presumably are related to basal topography (or slipperiness variations). Or maybe the assumption is that observed ice velocities reflect such … but is this correct?

AC: Yes, our assumption is that ice surface velocities observed from space reflect the variation of stress over the study area. This is probably reasonable at coarse temporal scales; the von-Mises stress field derived from the MEASURES monthly ice-flow velocity maps we use in this study clearly show larger portions of study area under high surface stress between July and August, which coincides with the maximum number of hydrofracture-induced SGL drainage events and rapid meltwater discharge events (through crevasses and moulins in the non-lake areas) occur. We agree that our approach is not strictly ‘correct’ at finer temporal scales – for example Christoffersen et al. (2020) show that basal lubrication due to the injection of surface runoff in the subglacial environment causes transient, accelerated ice-flow velocities, and thus more crevasses and lake drainage in the downstream areas. Again, this is something that we intend to develop in the future, once we have coupled with an ice sheet model and are thus able to force our hydrology model with velocity/stress fields at finer timescales. We have text to the revised manuscript to reflect this.

RC 1: In the quasistatic fracture propagation model described by Eqn 5, does conservation of water determine b?

AC: In Eqn 5, at every timestep (i.e., daily) we update the value of ‘b’ due to newly accumulated meltwater runoff. This is done by adjusting the total volume of accumulated surface meltwater runoff (at the corresponding crevassed cell) to the dimensions of the crevasse, and then updating the value of d (i.e., the crevasse depth) within a loop (i.e., by solving Eqn, 5) till K_I<K_IC. During our model runs, we found out that in every timestep, the new meltwater runoff that has accumulated at a crevassed cell will drive the crevasse deep enough for all of it to get squeezed inside the crevasse. This is similar to the findings of Krawczyinski et al. (2009), on whose work we based our model. They concluded that the crevasse depth constantly grows such that the level of the water inside the crevasse (i.e., b in our model runs) always stays below the crevasse’s surface.

RC 1: “Roughness” is poorly defined. Why 0.25? What is its range? I note that if f_R = 0, the assumed channel incision is zero and also lake height never decreases. This feels unphysical (melting can still happen even if water isn’t flowing if ice is at melting point). Please derive/justify this parameter.

AC: Roughness is poorly defined in the literature and developing a new definition is beyond the scope of this study. Merlind et al. (2006) express roughness as a function of the Manning roughness coefficient (n) and the hydraulic radius (R):

f_R=8gn²R^-1/3

For Greenlandic supraglacial streams, Merlind et al. (2006) estimated the value of n to be in the range 0.036 and 0.058 m^-1/3 with an average value of 0.050 m^-1/3, giving f_R values in the range of 0.2-1.5.  However, Kingslake et al. (2015) empirically derived 0.25 as the value of f_R for simulating channel incision for an SGL located in Greenland whose dimensions were acquired in situ by Giorgiou et al. (2009) and so we chose to use this value in keeping with their study. We checked the sensitivity of our model to Roughness by varying f_R between 0.25 and 1. The upper bound of 1 was derived by substituting the maximum value of daily MAR runoff (per 100 m-by-100 m DEM cell) for our study area which was ~200 mm, in the place of R in the above-mentioned equation. We found that the modelled results were not significantly different.

RC 1: Section 4.4: Does channelized or unchannelized water ever refreeze in this model? Isn’t there abundant evidence for refreezing of meltwater in firn?

AC: We assume that water flow in channels is transient, and that after meltwater production ceases, the channels become dry. Outside of channels, by using RCM estimates of runoff rather than melt, we assume that ‘refreezing of meltwater in firn’ is already taken care of in the MAR model.

RC 1: Paragraph 440: Is there data that validates the result of 100m-10s of km transport in a day? This translates to velocity of ~1 mm/s for 100 m/day, versus to ~0.1 m/s for 10km/day. There are numerous studies that have measured supraglacial stream discharge that you could compare with, for instance Smith et al., (2015).

AC: In our model, the water routing module transports water over distances of 100 m - 40 km day^-1. This amounts to a velocity range of 0.001 ms^-1- 0.462 ms^-1 and it overlaps with the range of 0.2-9.4 ms^-1, that was determined by Smith et al. (2015) from WorldView imagery. However, we note that the velocities were estimated in situ by Smith et al. (2015) and our estimates are averaged over a spatial resolution of 100 m.

We have added text to the revised manuscript giving our average flow speeds and comparing them to Yang et al. (2018) and Smith et al. (2015). We have also added Yang et al. (2018) and Smith et al. (2015) to the reference list.

References:

Smith, L. C., Chu, V. W., Yang, K., Gleason, C. J., Pitcher, L. H., Rennermalm, A. K., Legleiter, C. J., Behar, A. E., Overstreet, B. T., Moustafa, S. E., Tedesco, M., Forster, R. R., LeWinter, A. L., Finnegan, D. C., Sheng, Y., and Balog, J.: Efficient meltwater drainage through supraglacial streams and rivers on the southwest Greenland ice sheet, Proc. Natl. Acad. Sci., 112, 1001-1006,25 https://doi.org/10.1073/pnas.1413024112, 2015

Yang, K., Smith, L. C., Karlstrom, L., Cooper, M. G., Tedesco, M., van As, D., ... & Li, M. (2018). A new surface meltwater routing model for use on the<? xmltex\newline?> Greenland Ice Sheet surface. The Cryosphere, 12(12), 3791-3811.

Georgiou, S., Shepherd, A., McMillan, M., and Nienow, P.: Seasonal evolution of supraglacial lake volume from ASTER imagery. Annals of Glaciology, 50(52), 95–100, doi: 10.3189/172756409789624328, 2009

Kingslake, J., NG, F., and Sole, A.: Modelling channelized surface drainage of supraglacial lakes. Journal of Glaciology, 61, 225, doi: 10.3189/2015JoG14J158, 2015.

van der Veen, C.J.: Fracture propagation as means of rapidly transferring surface meltwater to the base of glaciers. Geophysics Research Letters, 34(1), L01501, doi: 10.1029/2006GL028385, 2008.

Colbeck, S. C.: The physical aspects of water flow through snow, in Advances in Hydroscience, vol. 11, edited by V. T. Chow, pp. 165–206, Academic, San Diego, Calif, 1978.

Mernild SH, Hasholt B and Liston GE (2006) Water flow through Mittivakkat Glacier, Ammassalik Island, SE Greenland. Geogr. Tidsskr., 106(1), 25–43 (doi: 10.1080/00167223.2006.10649543)

Karlstrom, L., & Yang, K. (2016). Fluvial supraglacial landscape evolution on the Greenland Ice Sheet. Geophysical Research Letters, 43(6), 2683-2692.

Crozier, J., Karlstrom, L., & Yang, K. (2018). Basal control of supraglacial meltwater catchments<? xmltex\break?> on the Greenland Ice Sheet. The Cryosphere, 12(10), 3383-3407.

Gudmundsson, G. H. (2003). Transmission of basal variability to a glacier surface. Journal of Geophysical Research: Solid Earth, 108(B5).

Krawczynski , M.J., Behn, M.D., Das, S.B., and Joughin, I.: Constraints on the lake volume required for hydro-fracture through ice sheets. Geophysical Research Letters, Vol. 36, L10501, doi:10.1029/2008GL036765, 2009.

We thank the Anonymous Referee #2 for their comments on our paper titled ‘A new model for supraglacial hydrology evolution and drainage for the Greenland ice sheet (SHED v1.0)’. The comments were constructive and have contributed in making this paper even better. In our response below, ‘RC 2’ stands for Comment by Referee #2 and ‘AC’ stands for Author Comments.

RC 2: Line 93: How are the hydrological catchments of regions 1, 2 and 3 calculated/where does the data about what the catchments are come from?

AC: For Region 1 and Region 2, the daily discharge data was acquired from the Kuussuup gauging station  and the Tasersiaq gauging station respectively, both of which are set up and maintained by the Asiaq Greenland Survey. This has been added in the caption to figure 1, in the revised manuscript.

Hydrological catchments for these regions are defined using the Arctic DEM (Porter et al., 2018) at a spatial resolution of 10 m using the Hydrology toolbox in ArcGIS. This was done by first filling sinks in the DEM, and then identifying all the cells that could potentially drain to the DEM cell at the location of each of the two gauging stations. The basin boundary is therefore the outermost ring of cells that meet this condition. The catchment was partitioned into ice and ice-free areas, by using the BedMachine V4 ice thickness map (Morlighem et al., 2017) of GrIS as a mask. This has now been stated in our revised manuscript.

For demonstrating the performance of the Lake refreezing module, we chose an arbitrarily large area (i.e., Region 3) surrounding the supraglacial lake (SGL) we use to test our model (Figure 1). The area was chosen to be large enough to capture all the influx of meltwater runoff into the SGL. Again, this is stated in the current version of our paper.

RC 2: Line 155: Here, the model is described as having 3 modules plus a separate model for lake refreezing. In figure 2, the model is described as having 4 modules, with lake refreezing being one of them. It would be good to have some consistency here, using the same number of modules in both places.

AC: Initially we had coded the refreezing module separately in Matlab. We have now coded this into Fortran as part of the main model, so there are now 4 module, and have rewritten the revised paper to read as follows:

‘Our model comprises four modules: 1) Supraglacial Routing and Lake Filling; 2) Hydrofracture; 3) Slow Lake Drainage and 4) Lake freezing module.’

RC 2: Line 164-166: I would appreciate some detail on how the PDC and the corresponding flow path are calculated.

AC: We have added a detailed description of this in Section 1 of the Appendices.

RC 2: Line 164-176: It sounds like quite a lot is being calculated at each timestep in the supraglacial routing module: the PDC and flow path for each DEM cell, the discharge and travel time between each pair of cells on the flow path. Would it be at all possible/useful to calculate all these things outright at the start (once, not every timestep), then apply it to the specific meltwater distribution at each timestep? This seems like it could be reasonable if the surface elevation isn’t changing very much over the model run (is this the case?), and could save a lot of computation. I expect that there might be a good reason for not doing this, but wondered whether it is something the authors have thought about.

AC: We agree with the reviewer that a lot of computation is performed at every timestep and that in this version of the model it would be more computationally efficient to perform these computations once only. However, we wanted to develop a hydrological model that can be coupled to a more complex ice-sheet model in the future. A coupled hydrology-ice-sheet model will update the DEM surface every model timestep. Therefore, we decided that the entire supraglacial routing module should be run every timestep. As a result, we estimate PDC, flowpath, discharge and travel time between each pair of cells along the flowpath, in every timestep.

RC 2: Line 172-174: What is the justification behind moving all the water from a DEM origin cell to the same new cell, rather than being distributed between multiple cells?

AC: In our paper, the hydrological model was designed by assuming that every SGL will have a unique surface catchment area (SCA) that feeds meltwater runoff into it. Consequently, the study area can be thought of as a combination of mutually exclusive SCAs (e.g., Arnold et al., 2010; Banwell et al., 2012; Koziol et al., 2017). In order to preserve this property of mutual exclusivity, we assumed that all the runoff water from a DEM origin cell will be transferred to a single destination cell only. This is stated in the revised version of the manuscript.

RC 2: Line 221-224: It is not clear to me how the crevasse depth d is being calculated here. In line 221, it says that ‘Equation 5 is solved for every crevassed cell...’, but what is being solved for: K_I or d? The way I understand it is that, at each timestep, equation 5 is used to calculate K_I, and if K_I>=K_IC then the crevasse propagates, with the new value of d in that timestep being calculated somehow, but I’m not sure how this new d is calculated. Is this correct? If so, how is d calculated? Or is it the case that, at each timestep, K_I is calculated using equation 5, then (i) if K_I<K_IC, d doesn’t change, and (ii) if K_I>=K_IC, the crevasse propagates with the new value of d for that timestep being set to whatever value of d gives K_I=K_IC? Some clarification on these points would be great.

AC: In every timestep, we estimate K_I from Eq. 5. If K_I>=K_IC i.e., 150 kPam^0.5, we increment value of d by 0.001 m (if there is sufficient water in the lake) and then re-estimate K_I from Eq. 5 using the updated value of d. This procedure is followed till either all of the water is contained in the crack and K_I gets less than K_IC or, d becomes equal to the DEM cell’s ice thickness. We have provided a better description of this process in our revised manuscript.

RC 2: Eqs. (6) and (9): I would appreciate a little intuition in the text as to where these equations come from, eg. why there is a power of 1.5.

AC: We agree with the reviewer and we have added a derivation of the expression in section 2 of the Appendices.

RC2: Line 271-272: It says here that H_c is determined both from lake geometry and by solving eqs. (6)-(9). Is this correct?

AC: Yes. The initial value of H_c is estimated from the lake geometry. This done by differencing the lake bottom elevation from lake outlet elevation for every SGL. Evolution of H_c is estimated by solving eqs. (6)-(9). This is now clarified more clearly in the revised manuscript.

RC 2: Line 272-273: Please clarify and explain where the overflow water from a SGL gets drained to. For example, when you say that the overflow water from the SGL is ‘routed as runoff via a meltwater channel downstream’, is this a separate process to module 1 (if so, what is the process?), or does the overflow water get immediately added to the supraglacial routing module?

AC: In a given timestep, the overflowing water is added to the destination cell as meltwater runoff and is then available for further routing through module 1 in the next timestep. This is now clarified more clearly in the revised manuscript.

RC 2: Line 279-281: Have you thought about whether this 1D approach of modelling only the deepest cell in the lake leads to under/over estimates of refreezing?

AC: Our estimates of lid-on and lid-break-up dates matched well with that observed from satellite imagery and so we think our approach is reasonable.

RC 2: Line 290: I have some concern here that some of the shortwave radiation may be double counted, being absorbed both at the surface and in the interior of the lake. The albedo \alpha is used to partition the shortwave radiation into a component \alpha*SWR that is reflected by the lake, and a component (1-\alpha)*SWR that is absorbed by the lake. Here, the surface energy balance (eq. (10)) features the term (1-\alpha)*SWR, suggesting that all of the shortwave radiation that the lake absorbs is absorbed at the surface. This is inconsistent with having a component of shortwave radiation FSW in eq. (14) that penetrates through the lake surface (being absorbed internally within the lake). To remove this inconsistency, I suggest including a parameter (eg. \chi in Hoffman et al 2014 (Journal of Glaciology), Law et al 2020 (Journal of Glaciology), Woods and Hewitt 2023 (The Cryosphere)) to account for the fact that not all of the non-reflected shortwave radiation is absorbed at the surface – some penetrates through the lake surface (as you mention in line 318). That is, if the total shortwave radiation absorbed by the ice is (1-\alpha)*SWR, then we can say that a proportion \chi of that is absorbed at the surface and the rest is absorbed below the lake surface (in the form of the penetrating radiation FSW), so the term in the SEB (10) would be \chi*(1-\alpha)*SWR and the below-surface absorbed radiation in the 1D eq. (14) would be FSW=(1-\chi)*(1-\alpha)*SWR. This avoids any shortwave radiation being double counted.

AC: We agree with the reviewers – it seems we made a typo in the paper. We have accordingly changed the formula in Eq. 10 and can confirm that the model code reflects the new formulation.

RC 2: Line 291-295: I would suggest listing the terms here in the same order than they appear in eq. (10).

AC: We agree with the suggestion and have implemented this suggestion.

RC 2: Line 315: Related to my line 290 comment, how do you calculate FSW?

AC: We have added more detail about this calculation into section 3 of the Appendix, as well as adding brief detail to the revised manuscript.

RC 2: Line 354: Why a value of 0.1m? Have other values been tried? Does this have much effect on the results?

AC: Previous studies such as Buzzard et al. (2018) have concluded that when a thin ice lid (< 0.1 m) forms on the SGL surface, it has a tendency to transiently disintegrate and then reform. This process continues till the air temperatures get low enough for the thin ice-lid to grow and become thick enough to resist disintegration. At this point, the ice-lid is classified as stable. Following Buzzard et al. (2018) and Law et al. )2020), we then define our lid-on date as the date when the lid thickness is positive.

We tested this assumption by also using 0.2 m as the threshold and the results did not differ much from the former. The lid-on date in 2016 moved back by a day but the lid-on and lid-break-up dates for other years remained the same. This is summarized more clearly in the revised manuscript.

RC 2: Line 370: It might be worth mentioning somewhere in section 4.4.1 that equation (14) is solved in the lake and (19) at the bottom, in addition to showing it in fig. 5.

AC: We agree with the suggestion and in the revised manuscript, we have added this sentence to the revised manuscript.

RC 2: Line 413-414: Why do these uncertainties affect the latter part of the melt season in particular?

AC: The uncertainties arise in the later part of the melt season due to:

1. Availability of more meltwater runoff as opposed to other months.
2. Presence of cloud over the study area leads to an under-estimation of satellite-imagery-measured daily lake area.

We have included this to the revised manuscript.

RC 2: Line 420: Why was a 5 m channel width chosen? Are there any observations to support this choice?

AC: We do not have any observations that can quantify the width of the supraglacial meltwater channels. Therefore, we chose the value of 5 m because this as found to be the most representative value of the range of channel widths observed in Western Greenland through WorldView imagery (Koziol et al., 2017). We did do a sensitivity analysis using 2 m as the channel width but the results were not very different from the experiment where, so the channel width was chosen as 5 m. This is explained in our paper

RC 2: Line 540: Is it consistent to define the modelled lid break-up date as the day when the ice lid thickness becomes zero, when a value of 0.1 m is used for defining the lid-on date? Also, is it consistent with the ‘30% exposed water surface’ definition for the satellite images? Related to this, I think it would be worth mentioning in section 4.4.1 what the condition for switching back to no lake ice is, when the relevant part of the model is first introduced.

AC:. In the revised manuscript, we have clarified that lid-on date as the date when the model first simulates 0.1 m of ice on the lake’s surface. We found out that with this assumption of ‘30% exposed water surface’ our results matched well with that observed.

RC 2: Figure 1: I suggest making clear what the axes mean (ie. longitude and latitude).

AC: We agree with the reviewer and now we have changed the figure. The new figure has labelled axes.

RC 2: Figure 2: Overall, this was very helpful for my understanding of the model – thank you for including it. One small suggestion I have is to put the second and third boxes in module 4 next to each other (horizontal) rather than vertically stacked, to make it clearer that it you only go through one of these two boxes in each cycle.

AC: Thank you for the suggestion, we have edited the figure accordingly.

RC 2: Figure 7: The elevation contours look solid, rather than dashed like the legend suggests.

AC: We agree with the reviewer and have amended the figure legend accordingly.

RC 2: Line 89: ‘were’ not ‘where’.

AC: Agreed. Changed

RC 2: Line 189: I assume ‘v^2’ is meant to say ‘\sigma_v^2’.

AC: Agreed. Changed

RC 2: Line 191: ‘\sigma_v^2’ not ‘\sigma_v 2’.

AC: Agreed. Changed

RC 2: Line 451: I think this should say ‘Figure 7’, not 8. Also, the year in the figure caption (2019) and in the text (2015) disagree.

AC: Agreed. Changed

RC 2: Line 461: I think this should say ‘Figure 7’, not 8.

AC: It is actually Figure 6c. Changed now.

RC 2: Line 467: I think this should say ‘Figures 8a and 8b’, not 10.

AC: Agreed. Changed

RC 2: Line 772: There’s a missing bracket at the start of ‘blue box’.

AC: Agreed. Changed

References:

Morlighem M. et al., (2017), BedMachine v3: Complete bed topography and ocean bathymetry mapping of Greenland from multi-beam echo sounding combined with mass conservation, Geophys. Res. Lett., 44, doi:10.1002/2017GL074954, http://onlinelibrary.wiley.com/doi/10.1002/2017GL074954/full

Porter, Cl., Morin, P., Howat, I., Noh, M., Bates, B., Peterman, K., Keesey, S., Schlenk, M., Gardiner, J., Tomko, K., Willis, M., Kelleher, C., Cloutier, M., Husby, E., Foga, S., Nakamura, H., Platson, M., Wethington, M,; Williamson, C., Bauer, G., Enos, J., Arnold, G., Kramer, W., Becker, P., Doshi, A., D'Souza, C., Cummens, P., Laurier, F., and Bojesen, M.: ArcticDEM, Harvard Dataverse. v2.0, doi:10.7910/DVN/ OHHUKH, 2018

Arnold, N. S.: A new approach for dealing with depressions in digital elevation models when calculating flow accumulation values, Prog. Phys. Geogr., 34(6), 781–809, doi:10.1177/0309133310384542, 2010

Banwell, A.F., Arnold, N., Willis, I., Tedesco, M., and Ahlstrom, A.: Modelling supraglacial water routing and lake filling on the Greenland Ice Sheet. J. of Geophysical Research. 117, doi:10.1029/2012JF002393, 2012.

Kingslake, J., NG, F., and Sole, A.: Modelling channelized surface drainage of supraglacial lakes. Journal of Glaciology, 61, 225, doi: 10.3189/2015JoG14J158, 2015.

Koziol, C., Arnold, N., Pope, A., and Colgan, W.:Quantifying supraglacial meltwater pathways in the Paakitsoq region, West Greenland. Journal of Glaciology, 63(239), 464-476, 2017.

We thank the Anonymous Referee #2 for their comments on our paper titled ‘A new model for supraglacial hydrology evolution and drainage for the Greenland ice sheet (SHED v1.0)’. The comments were constructive and have contributed in making this paper even better. In our response below, ‘RC 2’ stands for Comment by Referee #2 and ‘AC’ stands for Author Comments.

RC 2: Line 93: How are the hydrological catchments of regions 1, 2 and 3 calculated/where does the data about what the catchments are come from?

AC: For Region 1 and Region 2, the daily discharge data was acquired from the Kuussuup gauging station  and the Tasersiaq gauging station respectively, both of which are set up and maintained by the Asiaq Greenland Survey. This has been added in the caption to figure 1, in the revised manuscript.

Hydrological catchments for these regions are defined using the Arctic DEM (Porter et al., 2018) at a spatial resolution of 10 m using the Hydrology toolbox in ArcGIS. This was done by first filling sinks in the DEM, and then identifying all the cells that could potentially drain to the DEM cell at the location of each of the two gauging stations. The basin boundary is therefore the outermost ring of cells that meet this condition. The catchment was partitioned into ice and ice-free areas, by using the BedMachine V4 ice thickness map (Morlighem et al., 2017) of GrIS as a mask. This has now been stated in our revised manuscript.

For demonstrating the performance of the Lake refreezing module, we chose an arbitrarily large area (i.e., Region 3) surrounding the supraglacial lake (SGL) we use to test our model (Figure 1). The area was chosen to be large enough to capture all the influx of meltwater runoff into the SGL. Again, this is stated in the current version of our paper.

RC 2: Line 155: Here, the model is described as having 3 modules plus a separate model for lake refreezing. In figure 2, the model is described as having 4 modules, with lake refreezing being one of them. It would be good to have some consistency here, using the same number of modules in both places.

AC: Initially we had coded the refreezing module separately in Matlab. We have now coded this into Fortran as part of the main model, so there are now 4 module, and have rewritten the revised paper to read as follows:

‘Our model comprises four modules: 1) Supraglacial Routing and Lake Filling; 2) Hydrofracture; 3) Slow Lake Drainage and 4) Lake freezing module.’

RC 2: Line 164-166: I would appreciate some detail on how the PDC and the corresponding flow path are calculated.

AC: We have added a detailed description of this in Section 1 of the Appendices.

RC 2: Line 164-176: It sounds like quite a lot is being calculated at each timestep in the supraglacial routing module: the PDC and flow path for each DEM cell, the discharge and travel time between each pair of cells on the flow path. Would it be at all possible/useful to calculate all these things outright at the start (once, not every timestep), then apply it to the specific meltwater distribution at each timestep? This seems like it could be reasonable if the surface elevation isn’t changing very much over the model run (is this the case?), and could save a lot of computation. I expect that there might be a good reason for not doing this, but wondered whether it is something the authors have thought about.

AC: We agree with the reviewer that a lot of computation is performed at every timestep and that in this version of the model it would be more computationally efficient to perform these computations once only. However, we wanted to develop a hydrological model that can be coupled to a more complex ice-sheet model in the future. A coupled hydrology-ice-sheet model will update the DEM surface every model timestep. Therefore, we decided that the entire supraglacial routing module should be run every timestep. As a result, we estimate PDC, flowpath, discharge and travel time between each pair of cells along the flowpath, in every timestep.

RC 2: Line 172-174: What is the justification behind moving all the water from a DEM origin cell to the same new cell, rather than being distributed between multiple cells?

AC: In our paper, the hydrological model was designed by assuming that every SGL will have a unique surface catchment area (SCA) that feeds meltwater runoff into it. Consequently, the study area can be thought of as a combination of mutually exclusive SCAs (e.g., Arnold et al., 2010; Banwell et al., 2012; Koziol et al., 2017). In order to preserve this property of mutual exclusivity, we assumed that all the runoff water from a DEM origin cell will be transferred to a single destination cell only. This is stated in the revised version of the manuscript.

RC 2: Line 221-224: It is not clear to me how the crevasse depth d is being calculated here. In line 221, it says that ‘Equation 5 is solved for every crevassed cell...’, but what is being solved for: K_I or d? The way I understand it is that, at each timestep, equation 5 is used to calculate K_I, and if K_I>=K_IC then the crevasse propagates, with the new value of d in that timestep being calculated somehow, but I’m not sure how this new d is calculated. Is this correct? If so, how is d calculated? Or is it the case that, at each timestep, K_I is calculated using equation 5, then (i) if K_I<K_IC, d doesn’t change, and (ii) if K_I>=K_IC, the crevasse propagates with the new value of d for that timestep being set to whatever value of d gives K_I=K_IC? Some clarification on these points would be great.

AC: In every timestep, we estimate K_I from Eq. 5. If K_I>=K_IC i.e., 150 kPam^0.5, we increment value of d by 0.001 m (if there is sufficient water in the lake) and then re-estimate K_I from Eq. 5 using the updated value of d. This procedure is followed till either all of the water is contained in the crack and K_I gets less than K_IC or, d becomes equal to the DEM cell’s ice thickness. We have provided a better description of this process in our revised manuscript.

RC 2: Eqs. (6) and (9): I would appreciate a little intuition in the text as to where these equations come from, eg. why there is a power of 1.5.

AC: We agree with the reviewer and we have added a derivation of the expression in section 2 of the Appendices.

RC2: Line 271-272: It says here that H_c is determined both from lake geometry and by solving eqs. (6)-(9). Is this correct?

AC: Yes. The initial value of H_c is estimated from the lake geometry. This done by differencing the lake bottom elevation from lake outlet elevation for every SGL. Evolution of H_c is estimated by solving eqs. (6)-(9). This is now clarified more clearly in the revised manuscript.

RC 2: Line 272-273: Please clarify and explain where the overflow water from a SGL gets drained to. For example, when you say that the overflow water from the SGL is ‘routed as runoff via a meltwater channel downstream’, is this a separate process to module 1 (if so, what is the process?), or does the overflow water get immediately added to the supraglacial routing module?

AC: In a given timestep, the overflowing water is added to the destination cell as meltwater runoff and is then available for further routing through module 1 in the next timestep. This is now clarified more clearly in the revised manuscript.

RC 2: Line 279-281: Have you thought about whether this 1D approach of modelling only the deepest cell in the lake leads to under/over estimates of refreezing?

AC: Our estimates of lid-on and lid-break-up dates matched well with that observed from satellite imagery and so we think our approach is reasonable.

RC 2: Line 290: I have some concern here that some of the shortwave radiation may be double counted, being absorbed both at the surface and in the interior of the lake. The albedo \alpha is used to partition the shortwave radiation into a component \alpha*SWR that is reflected by the lake, and a component (1-\alpha)*SWR that is absorbed by the lake. Here, the surface energy balance (eq. (10)) features the term (1-\alpha)*SWR, suggesting that all of the shortwave radiation that the lake absorbs is absorbed at the surface. This is inconsistent with having a component of shortwave radiation FSW in eq. (14) that penetrates through the lake surface (being absorbed internally within the lake). To remove this inconsistency, I suggest including a parameter (eg. \chi in Hoffman et al 2014 (Journal of Glaciology), Law et al 2020 (Journal of Glaciology), Woods and Hewitt 2023 (The Cryosphere)) to account for the fact that not all of the non-reflected shortwave radiation is absorbed at the surface – some penetrates through the lake surface (as you mention in line 318). That is, if the total shortwave radiation absorbed by the ice is (1-\alpha)*SWR, then we can say that a proportion \chi of that is absorbed at the surface and the rest is absorbed below the lake surface (in the form of the penetrating radiation FSW), so the term in the SEB (10) would be \chi*(1-\alpha)*SWR and the below-surface absorbed radiation in the 1D eq. (14) would be FSW=(1-\chi)*(1-\alpha)*SWR. This avoids any shortwave radiation being double counted.

AC: We agree with the reviewers – it seems we made a typo in the paper. We have accordingly changed the formula in Eq. 10 and can confirm that the model code reflects the new formulation.

RC 2: Line 291-295: I would suggest listing the terms here in the same order than they appear in eq. (10).

AC: We agree with the suggestion and have implemented this suggestion.

RC 2: Line 315: Related to my line 290 comment, how do you calculate FSW?

AC: We have added more detail about this calculation into section 3 of the Appendix, as well as adding brief detail to the revised manuscript.

RC 2: Line 354: Why a value of 0.1m? Have other values been tried? Does this have much effect on the results?

AC: Previous studies such as Buzzard et al. (2018) have concluded that when a thin ice lid (< 0.1 m) forms on the SGL surface, it has a tendency to transiently disintegrate and then reform. This process continues till the air temperatures get low enough for the thin ice-lid to grow and become thick enough to resist disintegration. At this point, the ice-lid is classified as stable. Following Buzzard et al. (2018) and Law et al. )2020), we then define our lid-on date as the date when the lid thickness is positive.

We tested this assumption by also using 0.2 m as the threshold and the results did not differ much from the former. The lid-on date in 2016 moved back by a day but the lid-on and lid-break-up dates for other years remained the same. This is summarized more clearly in the revised manuscript.

RC 2: Line 370: It might be worth mentioning somewhere in section 4.4.1 that equation (14) is solved in the lake and (19) at the bottom, in addition to showing it in fig. 5.

AC: We agree with the suggestion and in the revised manuscript, we have added this sentence to the revised manuscript.

RC 2: Line 413-414: Why do these uncertainties affect the latter part of the melt season in particular?

AC: The uncertainties arise in the later part of the melt season due to:

1. Availability of more meltwater runoff as opposed to other months.
2. Presence of cloud over the study area leads to an under-estimation of satellite-imagery-measured daily lake area.

We have included this to the revised manuscript.

RC 2: Line 420: Why was a 5 m channel width chosen? Are there any observations to support this choice?

AC: We do not have any observations that can quantify the width of the supraglacial meltwater channels. Therefore, we chose the value of 5 m because this as found to be the most representative value of the range of channel widths observed in Western Greenland through WorldView imagery (Koziol et al., 2017). We did do a sensitivity analysis using 2 m as the channel width but the results were not very different from the experiment where, so the channel width was chosen as 5 m. This is explained in our paper

RC 2: Line 540: Is it consistent to define the modelled lid break-up date as the day when the ice lid thickness becomes zero, when a value of 0.1 m is used for defining the lid-on date? Also, is it consistent with the ‘30% exposed water surface’ definition for the satellite images? Related to this, I think it would be worth mentioning in section 4.4.1 what the condition for switching back to no lake ice is, when the relevant part of the model is first introduced.

AC:. In the revised manuscript, we have clarified that lid-on date as the date when the model first simulates 0.1 m of ice on the lake’s surface. We found out that with this assumption of ‘30% exposed water surface’ our results matched well with that observed.

RC 2: Figure 1: I suggest making clear what the axes mean (ie. longitude and latitude).

AC: We agree with the reviewer and now we have changed the figure. The new figure has labelled axes.

RC 2: Figure 2: Overall, this was very helpful for my understanding of the model – thank you for including it. One small suggestion I have is to put the second and third boxes in module 4 next to each other (horizontal) rather than vertically stacked, to make it clearer that it you only go through one of these two boxes in each cycle.

AC: Thank you for the suggestion, we have edited the figure accordingly.

RC 2: Figure 7: The elevation contours look solid, rather than dashed like the legend suggests.

AC: We agree with the reviewer and have amended the figure legend accordingly.

RC 2: Line 89: ‘were’ not ‘where’.

AC: Agreed. Changed

RC 2: Line 189: I assume ‘v^2’ is meant to say ‘\sigma_v^2’.

AC: Agreed. Changed

RC 2: Line 191: ‘\sigma_v^2’ not ‘\sigma_v 2’.

AC: Agreed. Changed

RC 2: Line 451: I think this should say ‘Figure 7’, not 8. Also, the year in the figure caption (2019) and in the text (2015) disagree.

AC: Agreed. Changed

RC 2: Line 461: I think this should say ‘Figure 7’, not 8.

AC: It is actually Figure 6c. Changed now.

RC 2: Line 467: I think this should say ‘Figures 8a and 8b’, not 10.

AC: Agreed. Changed

RC 2: Line 772: There’s a missing bracket at the start of ‘blue box’.

AC: Agreed. Changed

References:

Morlighem M. et al., (2017), BedMachine v3: Complete bed topography and ocean bathymetry mapping of Greenland from multi-beam echo sounding combined with mass conservation, Geophys. Res. Lett., 44, doi:10.1002/2017GL074954, http://onlinelibrary.wiley.com/doi/10.1002/2017GL074954/full

Porter, Cl., Morin, P., Howat, I., Noh, M., Bates, B., Peterman, K., Keesey, S., Schlenk, M., Gardiner, J., Tomko, K., Willis, M., Kelleher, C., Cloutier, M., Husby, E., Foga, S., Nakamura, H., Platson, M., Wethington, M,; Williamson, C., Bauer, G., Enos, J., Arnold, G., Kramer, W., Becker, P., Doshi, A., D'Souza, C., Cummens, P., Laurier, F., and Bojesen, M.: ArcticDEM, Harvard Dataverse. v2.0, doi:10.7910/DVN/ OHHUKH, 2018

Arnold, N. S.: A new approach for dealing with depressions in digital elevation models when calculating flow accumulation values, Prog. Phys. Geogr., 34(6), 781–809, doi:10.1177/0309133310384542, 2010

Banwell, A.F., Arnold, N., Willis, I., Tedesco, M., and Ahlstrom, A.: Modelling supraglacial water routing and lake filling on the Greenland Ice Sheet. J. of Geophysical Research. 117, doi:10.1029/2012JF002393, 2012.

Kingslake, J., NG, F., and Sole, A.: Modelling channelized surface drainage of supraglacial lakes. Journal of Glaciology, 61, 225, doi: 10.3189/2015JoG14J158, 2015.

Koziol, C., Arnold, N., Pope, A., and Colgan, W.:Quantifying supraglacial meltwater pathways in the Paakitsoq region, West Greenland. Journal of Glaciology, 63(239), 464-476, 2017.

We thank Referee #1 for their comments on our paper titled ‘A new model for supraglacial hydrology evolution and drainage for the Greenland ice sheet (SHED v1.0)’. The comments were constructive and have contributed in making this paper better. In our response below, ‘RC 1’ stands for Comment by Referee #1 and ‘AC’ stands for Author Comments.

RC 1: There are several ways to perform flow routing on a DEM (e.g., D8, Din, multislope etc). Which are you using and how do results vary if you change these? How do you account for DEM uncertainty?

AC: In our analyses, meltwater runoff was routed between the DEM cells, based on the principle of lowest lying neighbour according to which runoff from one DEM cell (i.e., origin) is transferred to a DEM cell that a) lies in the 3-by-3 neighbourhood of the origin and b) has the lowest surface elevation in the neighbourhood.

We also tested flow routing using the standard D8 algorithm according to which, runoff was transferred from the origin to a DEM cell that a) lies in the 3-by-3 neighbourhood of the origin and b) has the highest surface slope with respect to the origin. The results were not substantially different from our approach, and an additional advantage of our approach is that it is slightly more computationally efficient.

We have added text to explain about how we estimate flowpaths between a DEM cell and its Potential Destination Cell (PDC) in section 1 of Appendices and referenced the appendices in the main body of the revised manuscript.

The average uncertainty in elevation associated with the ArcticDEM product that we use in this study is ~0.1 m (Candela et al., 2017). Due to a paucity of in-situ elevation data for our region, we assume 0.1 m to be the vertical accuracy for the DEM of our study area.

We have added a sentence to the revised manuscript giving this uncertainty and have added the citation  Candela et al. (2017) to the reference list.

RC 1: Do you predict evolution of the underlying topography itself? Catchment boundaries, and supraglacial stream networks, likely evolve through the melt season. At the wavelengths of catchments, ice topography may be largely dictated by underlying bedrock topography (e.g., Gudmundsson , connected to stream networks explicitly by Crozier et al 2018). Thermal erosion in streams is rapid and thus dynamic reorganization of the supraglacial network (e.g., stream capture events and other significant/sudden changes to upstream drainage area at a point) are possible (Karlstrom and Yang, 2016).

AC: We simulate the evolution of the underlying surface topography in our modelled channels,by updating the surface elevation of the DEM cells (at every time step) as they are eroded in the ‘Lake overtopping and drainage module’ (i.e., Module 3). Thus it is possible for our model to simulate reorganization of the supraglacial meltwater channels when lakes drain laterally (e.g., Kalstrom and Yang, 2016). However, we do not account for change in surface elevation due to surface melting and ice-flow because we are not currently coupled with an ice flow model, but this would naturally be something we will investigate once we have a coupled model (which is a future goal beyond the scope of the current study). We have now added two sections of text to the revised manuscript to clarify this. This includes an explantion about how we simulate the re-organisation of meltwater channels as a result of rapid lateral lake drainage, as well as an explanation about our plan to move towards coupling with an ice-flow model in the future, so that we can simulate other observed processes (such as re-organisation of surface meltwater channels as a result of rapid vertical lake drainage events). We have also added Karlstrom and Yang (2016) to the reference list.

RC 1: Unchannelized flow through snow and firn is driven by hydraulic head above an impermeable substrate is an unconfined aquifer. If channels are known, one can calculate transit times based on the Dupuit approximation (e.g., Yang et al., 2018). It seems that porous flow through porous ice (not just snow) could be better accounted for – this represents a much slower transport than sheet flow implied by Manning. I was a bit surprised that you did not compare with the model of Yang et al., actually.

AC: When MAR predicts snow cover in a grid cell in our model, flow velocity and travel time is simulated as per Colbeck (1973). The travel time consists of two components. In the first, we estimate the time meltwater runoff takes to trickle down to the base of the snowpack, and in the second we estimate the time taken by the runoff to flow between DEM cells, at the snowpack’s base i.e., assuming saturated flow over impermeable substrate (by using Manning’s equation) (e.g., Arnold et al., 2020). We appreciate that weathered ice may also be porous, and that we do not account for this. However, given that the properties and extent of weathered ice are not known, we think that this is a reasonable approach, which has also been followed by previous studies (e.g., Koziol et al., 2017).

We note that our estimates of modelled flow velocity of meltwater runoff in supraglacial meltwater channels is in the range of 0.001 ms^-1- 0.462 ms^-1, in keeping with values suggested by Yang et al., 2018 (0.2-9.4 ms^-1). The reason our estimates are at the lower end of this range is likely because our estimates are modelled at a spatial resolution of 100 m, and those of Smith et al. (2015), which are used in Yang et al. (2018) are point measurements. We have added text to the revised manuscript stating our average flow speeds in comparison to Yang et al. (2018). We have also added Yang et al. (2018) to the reference list.

RC 1: Does variation of stress matter for assessing hydrofracture initation/propagation? This approach would seem to ignore contributions from bedrock/surface topography? Many depressions where lakes form presumably are related to basal topography (or slipperiness variations). Or maybe the assumption is that observed ice velocities reflect such … but is this correct?

AC: Yes, our assumption is that ice surface velocities observed from space reflect the variation of stress over the study area. This is probably reasonable at coarse temporal scales; the von-Mises stress field derived from the MEASURES monthly ice-flow velocity maps we use in this study clearly show larger portions of study area under high surface stress between July and August, which coincides with the maximum number of hydrofracture-induced SGL drainage events and rapid meltwater discharge events (through crevasses and moulins in the non-lake areas) occur. We agree that our approach is not strictly ‘correct’ at finer temporal scales – for example Christoffersen et al. (2020) show that basal lubrication due to the injection of surface runoff in the subglacial environment causes transient, accelerated ice-flow velocities, and thus more crevasses and lake drainage in the downstream areas. Again, this is something that we intend to develop in the future, once we have coupled with an ice sheet model and are thus able to force our hydrology model with velocity/stress fields at finer timescales. We have text to the revised manuscript to reflect this.

RC 1: In the quasistatic fracture propagation model described by Eqn 5, does conservation of water determine b?

AC: In Eqn 5, at every timestep (i.e., daily) we update the value of ‘b’ due to newly accumulated meltwater runoff. This is done by adjusting the total volume of accumulated surface meltwater runoff (at the corresponding crevassed cell) to the dimensions of the crevasse, and then updating the value of d (i.e., the crevasse depth) within a loop (i.e., by solving Eqn, 5) till K_I<K_IC. During our model runs, we found out that in every timestep, the new meltwater runoff that has accumulated at a crevassed cell will drive the crevasse deep enough for all of it to get squeezed inside the crevasse. This is similar to the findings of Krawczyinski et al. (2009), on whose work we based our model. They concluded that the crevasse depth constantly grows such that the level of the water inside the crevasse (i.e., b in our model runs) always stays below the crevasse’s surface.

RC 1: “Roughness” is poorly defined. Why 0.25? What is its range? I note that if f_R = 0, the assumed channel incision is zero and also lake height never decreases. This feels unphysical (melting can still happen even if water isn’t flowing if ice is at melting point). Please derive/justify this parameter.

AC: Roughness is poorly defined in the literature and developing a new definition is beyond the scope of this study. Merlind et al. (2006) express roughness as a function of the Manning roughness coefficient (n) and the hydraulic radius (R):

f_R=8gn²R^-1/3

For Greenlandic supraglacial streams, Merlind et al. (2006) estimated the value of n to be in the range 0.036 and 0.058 m^-1/3 with an average value of 0.050 m^-1/3, giving f_R values in the range of 0.2-1.5.  However, Kingslake et al. (2015) empirically derived 0.25 as the value of f_R for simulating channel incision for an SGL located in Greenland whose dimensions were acquired in situ by Giorgiou et al. (2009) and so we chose to use this value in keeping with their study. We checked the sensitivity of our model to Roughness by varying f_R between 0.25 and 1. The upper bound of 1 was derived by substituting the maximum value of daily MAR runoff (per 100 m-by-100 m DEM cell) for our study area which was ~200 mm, in the place of R in the above-mentioned equation. We found that the modelled results were not significantly different.

RC 1: Section 4.4: Does channelized or unchannelized water ever refreeze in this model? Isn’t there abundant evidence for refreezing of meltwater in firn?

AC: We assume that water flow in channels is transient, and that after meltwater production ceases, the channels become dry. Outside of channels, by using RCM estimates of runoff rather than melt, we assume that ‘refreezing of meltwater in firn’ is already taken care of in the MAR model.

RC 1: Paragraph 440: Is there data that validates the result of 100m-10s of km transport in a day? This translates to velocity of ~1 mm/s for 100 m/day, versus to ~0.1 m/s for 10km/day. There are numerous studies that have measured supraglacial stream discharge that you could compare with, for instance Smith et al., (2015).

AC: In our model, the water routing module transports water over distances of 100 m - 40 km day^-1. This amounts to a velocity range of 0.001 ms^-1- 0.462 ms^-1 and it overlaps with the range of 0.2-9.4 ms^-1, that was determined by Smith et al. (2015) from WorldView imagery. However, we note that the velocities were estimated in situ by Smith et al. (2015) and our estimates are averaged over a spatial resolution of 100 m.

We have added text to the revised manuscript giving our average flow speeds and comparing them to Yang et al. (2018) and Smith et al. (2015). We have also added Yang et al. (2018) and Smith et al. (2015) to the reference list.

References:

Smith, L. C., Chu, V. W., Yang, K., Gleason, C. J., Pitcher, L. H., Rennermalm, A. K., Legleiter, C. J., Behar, A. E., Overstreet, B. T., Moustafa, S. E., Tedesco, M., Forster, R. R., LeWinter, A. L., Finnegan, D. C., Sheng, Y., and Balog, J.: Efficient meltwater drainage through supraglacial streams and rivers on the southwest Greenland ice sheet, Proc. Natl. Acad. Sci., 112, 1001-1006,25 https://doi.org/10.1073/pnas.1413024112, 2015

Yang, K., Smith, L. C., Karlstrom, L., Cooper, M. G., Tedesco, M., van As, D., ... & Li, M. (2018). A new surface meltwater routing model for use on the<? xmltex\newline?> Greenland Ice Sheet surface. The Cryosphere, 12(12), 3791-3811.

Georgiou, S., Shepherd, A., McMillan, M., and Nienow, P.: Seasonal evolution of supraglacial lake volume from ASTER imagery. Annals of Glaciology, 50(52), 95–100, doi: 10.3189/172756409789624328, 2009

Kingslake, J., NG, F., and Sole, A.: Modelling channelized surface drainage of supraglacial lakes. Journal of Glaciology, 61, 225, doi: 10.3189/2015JoG14J158, 2015.

van der Veen, C.J.: Fracture propagation as means of rapidly transferring surface meltwater to the base of glaciers. Geophysics Research Letters, 34(1), L01501, doi: 10.1029/2006GL028385, 2008.

Colbeck, S. C.: The physical aspects of water flow through snow, in Advances in Hydroscience, vol. 11, edited by V. T. Chow, pp. 165–206, Academic, San Diego, Calif, 1978.

Mernild SH, Hasholt B and Liston GE (2006) Water flow through Mittivakkat Glacier, Ammassalik Island, SE Greenland. Geogr. Tidsskr., 106(1), 25–43 (doi: 10.1080/00167223.2006.10649543)

Karlstrom, L., & Yang, K. (2016). Fluvial supraglacial landscape evolution on the Greenland Ice Sheet. Geophysical Research Letters, 43(6), 2683-2692.

Crozier, J., Karlstrom, L., & Yang, K. (2018). Basal control of supraglacial meltwater catchments<? xmltex\break?> on the Greenland Ice Sheet. The Cryosphere, 12(10), 3383-3407.

Gudmundsson, G. H. (2003). Transmission of basal variability to a glacier surface. Journal of Geophysical Research: Solid Earth, 108(B5).

Krawczynski , M.J., Behn, M.D., Das, S.B., and Joughin, I.: Constraints on the lake volume required for hydro-fracture through ice sheets. Geophysical Research Letters, Vol. 36, L10501, doi:10.1029/2008GL036765, 2009.