Basal freeze-on generates complex ice-sheet stratigraphy

Large, plume-like internal ice-layer-structures have been observed in radar images from both Antarctica and Greenland, rising from the ice-sheet base to up to half of the ice thickness. Their origins are not yet understood. Here, we simulate their genesis by basal freeze-on using numerical ice-flow modelling and analyse the transient evolution of the emerging ice-plume and the surrounding ice-layer structure as a function of both freeze-on rate and ice flux. We find good agreement between radar observations, modelled ice-plume geometry and internal layer structure, and further show that plume height relates primarily to ice-flux and only secondarily to freeze-on. An in-depth analysis, performed for Northern Greenland of observed spatial plume distribution related to ice flow, basal topography and water availability supports our findings regarding ice flux and suggests freeze-on is controlled by ascending subglacial water flow. Our results imply that widespread basal freeze-on strongly affects ice stratigraphy and consequently ice-core interpretations.

depression of the melting point (PiDMP). Their major conclusion is that the features may be comprised of freeze-on and that the ice flux is the primary control on their height. The major claims of this paper are novel as several groups have been trying to isolate the source of these structures in the ice sheet. This work will be of interest to the ice rheology, ice core and eventually the ice sheet modeling communities. I have noted my concerns with the paper framing and the figures. I have concerns about Figure 1 have some suggestions on how to improve it. They should address why their hydrologic pathways differ from other published pathways. Others have indicated these features are along pathways ---it curious that theirs differ. "Thus, our spatial analysis supports the hypothesis of basal freeze-on due to ascending water for Greenland." In the summary section the authors note the presence of plumes in the regions of fast flow but do not address the distinct difference between Petermann and NEGIS. In Petermann there are plumes in the onset region while in NEGIS there are not. (line 140).

Other Comments
The authors omit mentions one of other primary difference between Antarctica and Greenlandthe well defined presence of water. In Antarctica there is a clear association with the basal plume and the subglacial water network. In Greenland the ability to map subglacial water is only recently been resolved. Currrently, we do not know the distribution of subglacial water in Greenland and its relationship to plumes. Freeze-on rates and subglacial water discharge The freeze-on rates necessary to create the plumes require a significant amount of subglacial meltwater. The authors justify the high melt rates by considering different heat fluxes and sources in the basal environment, and then calculate the subglacial water routes from the hydropotential. This seems like a simplified approach considering that several hydrological models have been developed for ice-flow models including Elmer/Ice . Even so, the simple approach does provide an estimate of freeze-on rates. It is not clear why this estimate of freezeon rates is not shown in Fig. 1 instead of the freeze-on index. Admittedly, the estimate is uncertain but by prescribing a high melt rate for the whole domain, the authors would get an upper boundary for the freeze-on rates. For the reader the freeze-on index is confusing and not easily translatable into freeze-on rate. The catchment area for the plume described and modelled in Fig. 2 is outlined in Fig. 1 as a large area close to the ice divide. In Fig. 1, however, only a few of the streamlines that initiate within the catchment area terminates or passes the plume. This suggests (to the reader) that the catchment area for the plume is in fact substantially smaller and by extension so is the amount of available water and thereby the freeze-on rates. Finally, studies have indicated that at least in the western part of North Greenland, the bed is likely to be frozen (MacGregor et al., 2016) thus it seems unlikely that large amounts of meltwater is passing through the subglacial system in this region. In contrast, in the eastern part the bed is likely thawed. This poses a problem for the freeze-on hypothesis since the plumes to the east of the ice divide would then be a result of a different process than the plumes west of the ice divide. Recommendation: Show (an upper bound for) freeze-on rates on the map, show that the outlined area is a catchment area for the plume in Fig. 2, and address the problem regarding the different basal thermal states of the ice sheet.
Ice-flow models The authors mention that the rheology of the plumes is different compared to meteoric ice. Unfortunately, this is not investigated further. The Elmer/Ice model does have the capability to resolve ice with different rheological properties and it would be extremely interesting to see the impact on ice-flow when a plume builds up with a different viscosity. Would the rheology contrast in itself be enough to sustain the plume (as claimed in other studies )? By including different ice rheologies the authors would be able to convincingly refute (or demonstrate) that rheology contrasts are important for the plume formation. Recommendation: Simulate plume formation and ice-flow invoking different ice rheologies. Fig. 1 is difficult to interpret. The contrast between orange, magenta and red (and blue) is not big enough. Overall, the figure is too busy. A legend would also help the reader. In Fig. 3d, it is not clear which line belongs to which axis (dotted line is h?). Line 23: assumed to be isochrones rather than "termed isochrones" Line 72: Distored = distorted? Line 109: Missing a space after full stop. Lines 121 and 122: The use of the >< symbol is not very helpful to the reader and could be phrased differently. Line 160: Presumably, "oldest ice" refers to the 1.5 million year old ice-core project but this should be stated more explicitly.

Minor comments
We would like to thank the Reviewers for their thorough and excellent comments and suggestions made.
We took the revision very seriously and believe that the findings in our manuscript are now stronger and clearer.
Notes to the following response: 'Lxx' refers to the Lines of the originally submitted manuscript (otherwise noted as 'new Lxx'), slanted black text is the reviewers comments, our comments are in orange, original text in black, and in blue is the changed or added text as found in the current manuscript.

Reviewers' comments:
Reviewer #1 (Remarks to the Author): The manuscript by Leysinger Vieli and al. present advances in the understanding of the formation of large plume-like features within ice layers, and previously identified from radar images both in Greenland and Antarctica. Although apparently common, the origin of these features remains unknown, and thus, their impact on ice dynamics and mass balance. This is a very interesting manuscript, with the most significant new result being that a plume height (several hundreds of meters) relates first to ice flux and only secondarily to freeze-on. This allows the authors to provide a physical explanation for large plume formations, even where basal freeze-on is low.
The manuscript is generally well written, but needs additional details in place. This would significantly improve clarity and thus help strengthen their results at times (see first comment below in particular).
Overall, this is great work that merits publication. My comments below mainly aim at improving clarity on the work done.
-L87-90: These lines introduce a key output from the study, yet there is not enough information for the new processes to be clear. I would suggest replacing "can be explained using mass-conservation principles and ice mechanics" with a paragraph providing more details on how these come into play.
This would also help clarify the link to the ratio R, which is not immediately obvious, and needs to be expanded upon (in other words, you show on Fig. 3c/d that plume height and ratio R are connected, but can you better explain why this is the case?).
Instead of referring only to the Methods we explain in more detail the principal of the flux relationship with regards to mass-conservation and ice mechanics. We continue from L88 (new L106) with the following: 1 In steady-state the total ice flux (Q) can be expressed by the balance flux (B), which is the sum of the accumulation rate (ȧ) and the freeze-on rate (ḟ ) integrated from the divide to the current location downstream, where for surface or basal meltȧ,ḟ < 0 (see Methods subsections 2.3 and 3, equations (2),(3) and (5)).
In summary, the total ice flux Q consists of ice flux arising from meteoric accumulation (M ) and from freeze-on (F ). In a vertical ice column in plane flow, the total ice flux Q is the sum of the meteoric ice flux M , which comes from the surface, and the freeze-on (or accreted) ice flux F , which originates at the bed (Fig. 5c). The relative height of each flux component is expressed by the partial flux through the respective ice column (see Methods section 3, equation (8)), and is a function of the relative contribution of M and B to the total ice flux Q and to the horizontal velocity distribution with depth (Methods equations (6) and (9)). (The relative contribution to total ice flux Q for each flux component is 1 = We define the fraction of the freeze-on flux component as ice-flux ratio which equals the normalised partial flux for an ice column of normalised height (h; see Methods section 3, equation (12)).
Along flow, the freeze-on flux F increases only over the accretion area for both the transient and the steady-state case, while the meteoric ice flux M is steadily increasing (Fig. 5c). This along-flow relationship between F and M leads, apart above accretion areas, to a decreasing flux ratio R (Fig. 5d). We find that the factor limiting the plume height is primarily given by the ice-flux ratio R.
-From Figure 1c / Supp Figure 1, it is difficult to see the good correspondence between the plume location and the places of crossing between subglacial water paths and elevated PHI, as it seems that quite a few do not fulfill both conditions. Would it be possible to specify in the text the proportion of plumes that are found to meet both criteria?
The analysis between Φ and the water path would indeed be nice if we knew where and at what rate the water is flowing. However, without this information this analysis adds not much information. This data set is complex and not very suitable for statistics, as plume information is available along radar transects only and all other information is calculated from surface and bed data, interpolated from RES data 11,12 ).
For freeze-on to occur we need the water flowing over a freeze-on area (elevated Φ) but the plume seen in the radar might be a bit further away along ice flow. Therefore, what we can do is to check if the ice at the plume comes from a nearby region crossing an accretion area. An argument that the observed plume relates to freeze-on, as the ice comes from an area with elevated Φ. We did that for the plumes ≥ 1/3H, < 1/3H and a random selection and added a new subsection 5.5 in the Methods section . The table is   added to the Supplementary Material (see also answer to Reviewer 3 on correlation between freeze-on   index and plumes).
-Overall, the method section is very long, which makes the reference to "Methods" not very helpful. Suggest adding numbering to the method subsections, and referring to specific method subsections throughout the manuscript.
We added numbering to the Methods subsections and now refer specifically to them throughout the text.
Furthermore, parts of the method section 'numerical models' have now been moved into the main text to improve clarity. However, the Methods section is still long.
-L 59/60: You should state here that you are using the SIA approximation, and why this is appropriate in your study. This could also be a good place to introduce the fact that you have additional runs to test the sensitivity of the model to the assumed simplified flow physics.
As mentioned above the sections removed from the numerical modelling section in Methods have been moved to that paragraph. We replaced L59/60 and added a sentence describing the model: We use a forward, three-dimensional, time-dependent, numerical ice-flow model conceived as a streamtube model (bed and surface are fixed in time; see Methods section 2.1) using finite differences. It uses simplified mechanics based on the shallow ice approximation (SIA), which ignores all horizontal stress gradients and is based on the assumption of small spatial stress gradients.
We replace L81-83 to continue with the explanation why the used models are appropriate and introduce We decided on adding the information directly into the text instead of providing a table. We added on L67 freeze-on area and rate (new L69): Following this, we simulate local freeze-on by adding mass using a freeze-on rateḟ = 0.8 m a −1 -of similar magnitude to the surface accumulation (ȧ = 0.1 m a −1 ) -over a basal accretion area 6 km along flow and 7 km across flow.
And we changed L84/85 adding the information on freeze-on rate and extent for Figure 3 to the following (new L100): In order to investigate the effect of ice flow on the plume shape, we apply an equal amount of basal freeze-on, using a freeze-on rateḟ = 0.42 m a −1 (surface accumulationȧ = 0.14 m a −1 ) along an extent of 5 km, at three different locations along the flow-line.
-This is a bit confusing: are you assuming in all your simulations that the ice moves purely by internal deformation? Or are you using eq. 6 from the methods section, which takes into account both basal sliding and internal deformation? If the former, this could be mentioned earlier on, e.g., as part of a model set-up description.
We add to L64 that we assume no basal sliding in all our calculations. On L94 where 'internal deformation' is introduced this is was already mentioned. L64 (new L64) now reads: The geometrical setting for a typical plume (Fig. 3a,b) is taken from observations (ice thickness, surface slope and distance to ice divide) 7,12,13 using simplified topography (uniformly inclined surface and a flat bed) while assuming no basal sliding.
-L. 110-111: You would expect that basal sliding increases for ice flow downstream. Shouldnt this also contribute to the smaller relative height of the plumes, according to Fig. 6? If that is the case, even if it is a secondary effect, you should include it in your interpretation. This is correct, relative plume height depends on the freeze-on flux F and the flux shape function (horizontal velocity distribution with depth). We add a sentence about the relation to the flux shape function to the sentence L110/111 (new L145) and refer to Figure 6.
Downstream, plume heights generally decrease due to the increasing meteoric ice flux M , which is a consequence of contributing accumulation area and flow convergence. In addition, basal sliding becomes more important downstream, affecting the distribution of horizontal velocity with depth, and thus further decreases the plume height (Fig. 6).
-L113-119: this paragraph seems a bit out of place. It needs better integration (to the following paragraph, presumably?).
Here the point is that we don't know much about the water at the base. However, if we look at PiDMP we know that the slope ratio is important and can amplify the freeze-on rate for any waterflux. We rewrote the paragraph (now two paragraphs) for better transition to the next paragraph, and also added a new subsection title.

Relation to basal freeze-on ratesḟ
In general, subglacial hydrology at the base of the ice sheet is poorly constrained as a consequence of the logistic difficulties of achieving sufficient observations. However, recent modelling studies suggest the existence of major meltwater pathways, subglacial lakes and, related to the mechanism described as the pressure-induced depression of the melting-point (PiDMP), widely occurring basal freeze-on 4,8 .
The predicted freeze-on rates are in general very low (< 1mm yr −1 ) and only locally larger rates (a few centimeters to meters) are obtained 4 . Substantial rates (tens of centimetres per year) have been previously suggested as a result of PiDMP, controlled by the adverse bed slope (α b ), surface slope (α s ) and the water flux (Q w ) 1 (see Methods section 4). The process of PiDMP leads to an increase in freezeon rateḟ both for increases in water flux Q w and bed slope α b , and to a decrease inḟ for an increase in surface slope α s 1,15,16 .
As we do not know the water flux Q w we use the relationship between the bed and surface slope to quantify the effect of PiDMP on the freeze-on rateḟ . Thus, we define the slope ratio S = α b /α s , and find that for slope ratios ranging between −2 ≥ S ≥ −11 the freeze-on rateḟ not only increases for increasing Q w but also for increasing α s ( Fig. 7 and Methods section 4.1, equation (15)). The effect of the slope ratios S on the freeze-on rateḟ , as shown in Figure 7a,b, ranges between a few centimeters (for both small α s and Q w ) to tens of centimeters (for large Q w ) and even meters (for both large α s and Q w ). It is readily conceivable that locally conditions arise, which lead to freeze-on rates comparable to the ones used in our experiments.
-L126-132: Quantifying the effects presented in this paragraph would strengthen its statement. Is the feedback mentioned important or not?
From the original Supplementary Figure 7 one can see how freeze-on rates change for a change in water flux for a given surface slope, or how the freeze-on rates change for a change in surface slope for a given water flux. We now changed the Figure so that it shows the specific slopes obtained from the full-system experiment. The changes in slope have a little effect on the freeze-on rate as it cancels each other out.
But more important is the effect on the hydraulic gradient. In the text we add the observed changes from the specific FS example. L126-130 now read (new L178): Our model experiment including full-flow physics and a freely evolving surface results in a reduced (increased) surface slope α s by 2 × 10 −4 upstream (downstream) of the freeze-on area and steadily increasing from 2.7 − 3.3 × 10 −3 over the freeze-on area (Fig. 4b,c), leading to a drop in slope ratio of nearly 2S along the bed slope α b , and hence leading to freeze-on ratesḟ that are slightly larger upstream than downstream assuming a constant water flux Q w (Fig. 7a,b, black line). The difference inḟ between the onset and end of the freeze-on area is a few percent and negligible over the whole area as they counterbalance each other. However, of significance is the effect of both spatial and temporal surface change 6 on the hydraulic gradient, possibly redirecting the water upstream or allowing water flow over steeper adverse slopes downstream. This complex feedback mechanism is expected to amplify and extend plume growth over initially unfavourable α b , while changing the hydraulic gradient. Furthermore, a distributed water system over an adverse slope α b steep enough for freeze-on to occur is likely to contribute to plume formation over a broader area 1 . Yes the effect is here also important. We changed the sentence on L141/142 and added a further sentence to the paragraph. It now reads (new L213): In fast flow areas, the large ice flux together with increased basal sliding (leading to a vertical distribution of the horizontal velocity resembling the plug flow case) prevent in general high rising plumes (Fig. 5).
An exception occurs for very large water fluxes and along-flow series of freeze-on areas, which leads to stacking of the advected plumes (Fig. 5b). Further along, in the ablation area, where surface melt reduces both the meteoric ice flux M and, as a consequence, the total ice flux Q (hence increasing R) and high availability of surface melt-water reaches the bed 2,10 , larger plumes are again possible.
We agree that the sentence is not entirely clear. We try to highlight how the different components, we looked at so far, influence each other and match the observeations. The sentence now reads (new L247): In summary, our arguments regarding the interplay between ice flux and PiDMP through the freeze-on rateḟ , which is controlled by the basal water flux Q w , the freeze-on index Φ, and the location of the observed plumes, support the hypothesis of basal freeze-on.
-Finally, there is a lack of discussion on how the new findings fit / relate to existing work on the formation of theses plumes. I would suggest adding details either in the intro and/or conclusion.
We added details by comparing our results to previous work throughout the result section and also added a new section on viscous contrast which produces structures that were discussed in recent literature. We put our results in context with newly published work. should be moved to that particular section.
We reorganised the original Figure 1  We regrouped the information for the panels to make it clearer and followed your suggestions. We took up your suggestion and changed caption and figure accordingly.
Reviewer #2 (Remarks to the Author): Review of Basal freeze-on at the origin of complex internal ice-layer stratigraphy.
Title: Might consider clearer title: Extent of basal freeze-on in Northern Greenland This suggested title focuses on the extent of basal freeze-on, however our focus is on the process of basal freeze-on. Therefore we believe the original title is expressing this aspect accurately.
Overview: This paper works to address systematically a very intriguing question -how much of the ice at the base of the Greenland Ice Sheet is refrozen basal water and how does the refreezing work. This is a nice integration of data analysis and modeling. The authors examine the pressure-induced depression of the melting point (PiDMP). Their major conclusion is that the features may be comprised of freeze-on and that the ice flux is the primary control on their height. The major claims of this paper are novel as several groups have been trying to isolate the source of these structures in the ice sheet. This work will be of interest to the ice rheology, ice core and eventually the ice sheet modeling communities. I have noted my concerns with the paper framing and the figures.
I have concerns about Figure 1 have some suggestions on how to improve it. They should address why their hydrologic pathways differ from other published pathways. Others have indicated these features are along pathways -it curious that theirs differ.
"Thus, our spatial analysis supports the hypothesis of basal freeze-on due to ascending water for Greenland." In the original paper we calculated the streamlines along the reversed hydraulic gradient in order to obtain the source region. This is why the hydrological pathway looks different than in other publications.
It was not stated clearly in the text. Now we plot the hydrological pathways in following the hydraulic gradient downwards, representing how the water flows from the seed point towards the ocean.  Near the ice divide the observed plumes are usually smaller than 1/10H, in the range of a few 100 meters. For these relative small plume heights the flux ratio R too is small when assuming flow purely by internal deformation. As the ice flux is low, only low water fluxes are needed to produce the required freeze-on rates (a few millimeters to centimeters). Since the surface slope is nearly flat, Φ is low too (< 10 −3 ), commonly leading to bed slopes that are too steep for water to flow uphill, and restricting plume formation. Away from the divide, the first large plumes (> 1/3H) start to form in regions where the basal water flux Q w is increasing, with the total ice flux Q and freeze-on index Φ still being generally small. Further along, towards the margin both Φ and Q w tend to increase (while Q is still moderate) and allow for larger plume heights. In regions of fast flow (large Q) the available basal water flux and Φ are in general too low to produce freeze-on ratesḟ that would sustain the large ice-flux ratio R required for large plumes. These observations fit well with the modelling results by Dow et al. 4 producing larger freeze-on rates along important water flow paths and towards the margin.

Other Comments
The authors omit mentions one of other primary difference between Antarctica and Greenland the well defined presence of water. In Antarctica there is a clear association with the basal plume and the subglacial water network. In Greenland the ability to map subglacial water is only recently been resolved.
Currrently, we do not know the distribution of subglacial water in Greenland and its relationship to plumes.
We agree that despite some effort (e.g. Livingston et al 2013; Dow et al, 2018) we do not know the distribution of subglacial water in Greenland. This is why we assume that we have water everywhere to map possible water paths for the case that there is water. But from ice core drilling we do know that there are regions at the base with water. So we know that there is water but we do not know where.
Figures: Figure 1 is terribly complicated and difficult to follow. I could recommend several revisions.  A plot of plume height versus freeze-on index is a nice idea. But it would not convey the relationship between freez-on index and plume height as the height depends on various factors. Ice flux, the flux of basal water, the slope ratio, the spatial area as well as the temporal extent of freeze-on influence the height of the plume. In areas of low ice flux large plumes can be obtained with relative little freezeon and therefore not requiring a large freeze-on index. However, towards the margin, with increasing ice flux, large freeze-on is required and therefore rather a larger freeze-on index. Only for the same condition a larger freeze-on index would lead to a higher plume. As the freeze-on index is obtained using information derived from bed and surface topography (slope, hydropotential) such a plot would not clearly highlight the relationship of freeze-on index and plume hight.
Sup. Figure  (see also answer to Reviewer #1). We find that if we seed within 1.5 km from the observed large plumes about 97% have a freeze-on area 5 km upstream, when calculated on a grid of 150 m. While for the small plumes 93% cover such an area. Doing the same with a random distribution of locations the fit is below 80%. Calculating the same on the larger grid of 1050 m as used in the paper figures we obtain similar results in the sense that the larger plumes give a higher match with freeze-on areas than for smaller plumes. This is most likley due to the fact that it is easier to distinct for large plumes if they rise from the bed or if they are advected. The plumes < 1/3H are also mapped when the RES layer structure of the meteoric ice shows the typical upward bending seen from the modelling experiment when active basal freeze-on is happening. We added the following  We include all plumes on the map, differencing clearer between large plumes and small plumes by choosing different symbols and sizes. Small plumes are not expected to behave differently, but as mentioned above the mapping is not as obvious as for larger ones. Our reason to leave out the small plumes was that we only had them consistently for 2010 − 2012. We now updated the mapping for all plumes for the years 2013 − 2014.

Freeze-on rates and subglacial water discharge
The freeze-on rates necessary to create the plumes require a significant amount of subglacial meltwater.
The authors justify the high melt rates by considering different heat fluxes and sources in the basal environment, and then calculate the subglacial water routes from the hydropotential. This seems like a simplified approach considering that several hydrological models have been developed for ice-flow models including Elmer/Ice (de Fleurian et al., 2014). Even so, the simple approach does provide an estimate of freeze-on rates.
We agree that our calculation of the subglacial water pathways, assuming that we have water everywhere, is a simplified approach. Using or even developing a hydrological model is beyond the scope of this paper. Especially as there is no subglacial data to validate the models. Therefore we prefer to know that our approach is simplified but gives us some indication where water would likely flow. A recent study  applied the subglacial hydrology model Glacier Drainage System (GlaDS) to assess the locations and rates of freeze-on by supercooling (freeze-on due to a pressure-induced depression of the melting-point (PiDMP)) for both distributed and efficient drainage networks. They found that freezeon occurred in many areas of the ice-sheet but obtained very low averaged freeze-on rates. However, locally freeze-on rates could reach substantial values (centimeters to meters) with the highest found in channels and the regions of higher freeze-on rates broadly matching the spatial distribution of our mapped plume-like features.
We added to the new paragraph starting on L135 (see comments to Reviewer #2) the following sentence: These observations fit well with the modelling results by Dow et al. 4 producing larger freeze-on rates along important water flow paths and towards the margin.
It is not clear why this estimate of freeze-on rates is not shown in Fig. 1 instead of the freeze-on index. Admittedly, the estimate is uncertain but by prescribing a high melt rate for the whole domain, the authors would get an upper boundary for the freeze-on rates. For the reader the freeze-on index is confusing and not easily translatable into freeze-on rate.
We chose by purpose to show the freeze-on index, as a simple multiplication with the water flux, ice density, gravitational acceleration and a division of the whole product by the volumetric latent heat of freezing (see Equation (19) and (20)) leads to the freeze-on rate. By assuming a constant basal melt rate we would not expect the water flux to be constant but rather expect to vary locally. Therefore the freeze-on index is of more general value and can be used with whatever water flux is appropriate for the region of interest.
The catchment area for the plume described and modelled in Fig. 2 is outlined in Fig. 1 as a large area close to the ice divide. In Fig. 1, however, only a few of the streamlines that initiate within the catchment area terminates or passes the plume. This suggests (to the reader) that the catchment area for the plume is in fact substantially smaller and by extension so is the amount of available water and thereby the freeze-on rates.
We adjusted the outlined catchment area above the plume by restricting it only to the streamlines that feed the plume. It was not clearly stated in the text that the streamlines were calculated to show where the water came from. The streamlines plotted together with our previous outline were calculated by streaming from the plume upwards along the hydrological gradient. Now we seed over the accretion length and 15 km upstream along ice flow using a 1050 m grid and calculating streamlines using steps of 5 km. The path of the streamlines is shown in the new overview Figure (Supplementary Fig. 2c). The previous figure used seeds at 30 km spacing with 3 km steps. As our method is not filling up sinks we need to step over them. The larger step size helps to cross regions where water is ponding and therefore reaches further upstream than a smaller step size. The outlined catchment area with a step of 5 km To obtain the source area we seed at 1050 m along the 6 km accretion length (observed in RES profile) and 15 km plume upstream along ice flow, and calculate streamlines along the reversed (upwards) hydraulic gradient using steps of 5 km ( Supplementary Fig. 2d).
Finally, studies have indicated that at least in the western part of North Greenland, the bed is likely to be frozen (MacGregor et al., 2016) thus it seems unlikely that large amounts of meltwater is passing through the subglacial system in this region. In contrast, in the eastern part the bed is likely thawed.
This poses a problem for the freeze-on hypothesis since the plumes to the east of the ice divide would then be a result of a different process than the plumes west of the ice divide.
We find that when using the outlines of MacGregor et al (2016)  We added a new overview figure (Supplementary Fig. 2b) showing streamlines seeded in the areas' likley thawed' and 'uncertain', which feed some of the plumes in the 'likely frozen' area with water.
Streamlines seeded only in the 'likely thawed' area ( Supplementary Fig. 2c) do leave a substantial amount of plumes in the Petermann and NEGIS region without water.
We added following in the methods section 5.1 (new L487): Comparing our mapped plumes with a composition of areas with a temperate bed 9 we find that most plumes are within areas of a 'likely thawed' or 'uncertain' bed, where water paths seeded from both areas reach most plumes ( Supplementary Fig. 2b,c).
Recommendation: Show (an upper bound for) freeze-on rates on the map, show that the outlined area is a catchment area for the plume in Fig. 2, and address the problem regarding the different basal thermal states of the ice sheet.
We addressed these points in our responses above.

Ice-flow models
The authors mention that the rheology of the plumes is different compared to meteoric ice. Unfortu- we defined the accreted ice with a contrasting viscosity to the meteoric ice. We performed two suites of experiments where we create a twofold and tenfold viscosity contrast between the meteoric and the accreted ice. The accreted ice is changed to being two or ten times stiffer (factor of 1/2 or 1/10) than the surrounding (meteoric) ice and changed to being two and ten times softer (factor of 2 or 10) than the meteoric ice. We find, as predicted by the kinematic theory, that over the accretion area the plume with stiffer freeze-on ice rises higher than the plume without viscosity contrast (since the ice can't shear easily it grows higher) and that the plume with softer ice rises less high (as it can shear easily it will move downstream). However, at a distance downstream of the accretion area both plumes with a viscosity contrast rise higher than for the plume without a contrast. The result of the softer ice is not entirely obvious -one way to explain it, is that while the soft ice shears easily it will be blocked by the stiffer ice downstream of the accretion area and has to rise over it. However, as it is a complex fluid dynamical problem it may therefore not have a simple explanation. What this experiment shows is that using a vis-cosity contrast leads to general higher plumes, therefore requiring smaller freeze-on rates (and therefore smaller water fluxes) to produce the observed plume-like structures. Further the results of the experiment suggests that from the plume shape we can gain information on the viscosity contrast.
We added a figure (Fig. 8) and following paragraph to the main text continuing on L144 (new L219):

Contrast in ice rheology
Plume-ice rheology is expected to differ from meteoric ice rheology 2 , reflecting different ice temperature and fabric, and consequently affecting the relationship between ice-flux ratio R and relative height h (equation (12)) for cases with internal deformation rather than plug flow. Large differences in rheological properties in near basal ice has been observed in Greenland 3 . A recent study by Wrona et al. 18 describing the geometry and morphology of near basal anomalous layer structures in RES profiles over the Gamburtsev mountains, East Antarctica, concludes that mechanical mixing between meteoric and accreted ice is triggered by rheology contrasts producing complex structures.
To investigate the effect of a rheology contrast between the meteoric and accreted ice on plume growth, we perform a simple experiment, using the FS model. We use the same experiment set up as described ice rises not as high as in the original case (Fig. 8a,d), which might be explained by the fact that it shears more easily. However, downstream of the accretion area the softer ice eventually rises higher than the original plume, which is simply advected downstream (Fig. 8b,e). A possible explanation could be that downstream of the accretion area the soft ice has to override the stiffer meteoric ice. Accretion of softer ice over a long time scale leads to instabilities along the advected plume margin (Fig. 8c,f) resembling the 'fingering' structures described by Wrona et al. 18 . Since a fluid infiltrating a higher viscosity fluid is a complex fluid dynamical problem, it might not have a simple explanation. Nevertheless, we can state that the plume height highly depends on the viscosity contrast, requiring smaller freeze-on rates to achieve the same plume height than stated above. We find that, whatever the rheological contrast (within reason), plume structures can be formed by the accretion process. Fig. 1 is difficult to interpret. The contrast between orange, magenta and red (and blue) is not big enough. Overall, the figure is too busy. A legend would also help the reader.

Minor comments
We tried to improve this figure by splitting it into several figures and keeping the originally Fig. 1c simpler with stronger colors using a coarser grid (1050 m instead of 600 m). See comments to Reviewer #1.
In Fig. 3d, it is not clear which line belongs to which axis (dotted line is h?).
The axis are valid for both lines -there is a different grid -so one can read both values along one line. Line 23: assumed to be isochrones rather than "termed isochrones" Changed this in text.
We corrected this.
Line 109: Missing a space after full stop.
Added space after full stop.
Lines 121 and 122: The use of the ¿¡ symbol is not very helpful to the reader and could be phrased differently.
We changed the sentence starting on L120 with (new L169): The dependence of PiDMP on S has the consequence that water can flow up steeper adverse bed slopes α b > 0 for a steeper inclined surface slope α s < 0 (or vice versa) , leading to a larger PiDMP effect owing to larger spatial gradients in overburden and thus melting point.
Line 160: Presumably, "oldest ice" refers to the 1.5 million year old ice-core project but this should be stated more explicitly.
We added to L159/160 with (new L266): Furthermore, plume growth alters the age-depth relationship, as bottom layers are lifted and compressed, while protected from basal melt, with important implications for ice-core drilling, interpretations, and the search for a site to core at least 1.5-million-year-old ice 14 .
Reviewers' comments: Reviewer #1 (Remarks to the Author): The authors have made thorough revisions to the manuscript, which I now find much clearer. Equally, the revisions to the figures are very helpful.
Overall, this is a very interesting manuscript, which will no doubt be of interest to the glaciology community.
I have just a few additional comments, which could be addressed by minor adjustment to the text.
-The first comment concerns the relationship between the plume location / elevated PHI / water paths: In their response, the authors explain that it is difficult to show the correlation between plume location and high PHI / water path, due to various uncertainties in the available datasets. Thus, I would suggest that this is reflected in the text by changing the word "reveals" on line 57 for something less assertive (e.g. Reviewer #2 (Remarks to the Author): Review of Revised Manuscript: Basal freeze-on at the origin of complex internal ice-layer stratigraphy.
In general the paper has been significantly improved in revision. The expanded text makes the work cleared and the figures are much better and easier to understand. My concerns about Figure  1 in particular have been addressed. The authors response to my questions about the difference in the plume distribution between Petermann and NEGIS catchments in northern Greenland clarifies the issue. Their statement that we assume that water is plentiful at the base of the Greenland ice sheet is adequate. The linkage between the main text and the supplementary figures is improved and the supplementary figures are improved. I would still suggest a clearer title such as: Contribution of basal freeze-on to complex ice sheet stratigraphy My argument would be that the "at the origin" is not the critical contribution of this paper. The authors have responded to the reviews comments adequately and I would suggest publishing the manuscript.
This is a very nice modelling study showing simulations of internal stratigraphy development under different conditions. However, I have some reservations about the conclusions of the paper. While I agree that the model results are similar to the observed layer stratigraphy I do not think that is a strong enough argument in itself. Especially since several other studies have successfully reproduced the layer stratigraphy invoking other processes. Thus, while the authors demonstrate that it is possible to grow the plume-like features with basal freeze-on, I am not convinced that it is the most probable process. Below, I outline my main concerns.
#1 Spatial correlation between plumes, freeze-on index and subglacial water paths Lines 51-53: "Analysing the spatial distribution of internal ice-layer structures with regard to Φ and water flow-path, reveals that large plumes are located along subglacial water paths crossing areas of elevated Φ ( Fig. 2 and Supplementary Fig. 2; see Methods section 5.1)." Lines 191-192: "A detailed spatial analysis for Greenland locates large plumes along possible water-paths crossing areas of high freeze-on indices…" Suppl. Table 1.
It is not obvious to me that there is a spatial correlation. From Fig. 2, areas of elevated freeze-on index occur in numerous areas. In some areas, they coincide with subglacial water paths and plumes, and in other areas they do not (see also attached figure). The attempt to correlate plume locations and freeze-on index in Suppl. Table 1 does not shed light on this. Considering that the freeze-on index is high in so many parts of the ice sheet, I would expect most streamlines to cross such an area at some point.
Lines 147-149: "In addition, basal sliding becomes more important downstream, affecting the distribution of horizontal velocity with depth, and thus further decreases the plume height (Fig.  6)." Lines 213-216: "In fast flow areas, the large ice flux together with increased basal sliding (leading to a vertical distribution of the horizontal velocity resembling the plug flow case) prevent in general high rising plumes (Fig. 5). An exception occurs for very large water fluxes and along-flow series of freeze-on areas, which leads to stacking of the advected plumes (Fig. 5b)." Observations from Petermann basin (e.g. Bell et al., 2014 and Panton and Karlsson, 2015) indicate that the plumes are increasing in size downstream. In Fig. 2, there is a large area with elevated freeze-on index but the plumes are downstream of this area. This seems to be in direct contrast with the modelling results (e.g. Fig. 5b where the plumes decrease in size downstream). Are the plumes observed in Petermann a result of "very large water fluxes"? And in this case, how large? Can the model reproduce the observations from Petermann?
Lines 199-202: "Away from the divide, the first large plumes (> 1/3H) start to form in regions where the basal water flux Qw is increasing, with the total ice flux Q and freeze-on index Φ still being generally small. Further along, towards the margin both Φ and Qw tend to increase (while Q is still moderate) and allow for larger plume heights." Based on the suppl. Fig. 2, I would argue that it is just as likely that the plumes form in areas where there is a transition from frozen to thawed conditions. #2 Amount of basal meltwater and freeze-on rates Line 154: Reference no. 27 and 28. Both references here are modelling studies and do not provide direct evidence of subglacial water or major pathways. Line 167-168: "It is readily conceivable that locally conditions arise, which lead to freeze-on rates comparable to the ones used in our experiments." Is it? This statement comes across as unsubstantiated. Is the slope ratio reasonable for the areas where the plumes are observed? Are the values for the water flux reasonable? This needs to be stated explicitly. Also, there is a difference between conditions such as these arising locally, and to them being present on a large scale across the entire ice sheet as evidenced by the extensive presence of the plumes. Lines 251-252: "…and broadly coincide with modelled freeze-on regions [28]" This statement is a bit misleading since the freeze-on rates in the study cited (Dow et al.) are at least an order of magnitude smaller than what is invoked here. Lines 262-263: "Therefore, statements about the mismatch between plume size and freeze-on rates [28] are not valid without considering as well the effect on the relationship between plume height and freeze-on rate both from the flux relationship and rheology contrast." I agree that looking at the freeze-on rate in itself is not adequate when examining plume heights. However, from the figures it seems that 10*softer ice roughly leads to a doubling in plume height.
Since the freeze-on rates are one to two orders of magnitude larger than the rates from Dow et al., using the freeze-on rates from the latter study would never lead to a plume of that height, regardless of whether or not ice-flux and rheology contrasts were included. It should be explicitly stated somewhere that the freeze-on rates from Dow et al., would only lead to X m (or x relative height) plumes. Otherwise the reader might think that the addition of ice-flux and rheology contrast in itself is enough.
#3 Frozen basal conditions Lines 194-197: "Near the ice divide the observed plumes are usually smaller than 1/10H, in the range of a few 100 meters. For these relative small plume heights the flux ratio R too is small when assuming flow purely by internal deformation. As the ice flux is low, only low water fluxes are needed to produce the required freeze-on rates (a few millimeters to centimeters)." Near the ice divide, radar observations indicate that the bed is likely frozen (as shown in Suppl. Fig. 2), thus the subglacial water flux must be zero and therefore also the freeze-on rates. It is a significant weakness in this manuscript that the presence of the plumes in the interior of the ice sheet cannot be explained by the freeze-on process but this weakness is not mentioned or discussed.

Response to the Reviewer's comments
Notes to the following response: 'Lxx' refers to the Lines of the originally submitted manuscript (otherwise noted as 'new Lxx', which refers to the manuscript with mark up), slanted black text is the reviewers comments, our comments are in red, original text in black, and in blue is the changed or added text as found in the revised manuscript.
We made a few general changes in this third revised version:  (v) Throughout the text we corrected typos and improved sentences.
The authors have made thorough revisions to the manuscript, which I now find much clearer. Equally, the revisions to the figures are very helpful.
Overall, this is a very interesting manuscript, which will no doubt be of interest to the glaciology community.
I have just a few additional comments, which could be addressed by minor adjustment to the text.
-The first comment concerns the relationship between the plume location / elevated PHI / water paths: In their response, the authors explain that it is difficult to show the correlation between plume location and high PHI / water path, due to various uncertainties in the available datasets. Thus, I would suggest that this is reflected in the text by changing the word "reveals" on line 57 for something less assertive -line 108: would suggest using "a uniform amount of basal freeze-on" Agreed -changed 'equal' to 'uniform'.
-line 109: the values of freeze rate (0.42 m/yr) and accumulation rate (0.14m/yr) seem very specific.

Can you explain better what drove this parameter choice?
We chose the freeze-on rate to be three times larger than the accumulation rate, as to produce large plumes comparable to observed ones that show larger effects in the deflection of the internal layers of the meteoric ice. It is clearer to highlight the relationship of the freeze-on rate to the surface accumulation instead of the exact values chosen. We replaced this in the text (new L132-135) In order to investigate the effect of ice flow on the plume shape, we apply an uniformequal amount of basal freeze-on, using a freeze-on rate three times larger than the surface accumulationḟ = 0.42 m a −1 (surface accumulationȧ = 0.14 m a −1 ) along an extent of 5 km, to produce distinctive plumes at three different locations along the flow-line.
We leave the definition of 'adverse slope' where it is but remove 'adverse' in the first sentence. It now reads (new L190-192): Substantial rates (tens of centimetres per year) have been previously suggested as a result of PiDMP, controlled by the adverse bed slope (α b ), surface slope (α s ), and the water flux (Q w ) when α b has opposite sign to α s 1 (see Methods section 4).

Reviewer #2 (Remarks to the Author):
Review of Revised Manuscript: Basal freeze-on at the origin of complex internal ice-layer stratigraphy.
In general the paper has been significantly improved in revision. The expanded text makes the work I would still suggest a clearer title such as: Contribution of basal freeze-on to complex ice sheet stratigraphy My argument would be that the at the origin is not the critical contribution of this paper. The authors have responded to the reviews comments adequately and I would suggest publishing the manuscript.
We see the point of the reviewer and agree that the suggested title is clearer and more accurate than our previous title. Prior to our study it was not expected that basal freeze-on can produce structures of the observed vertical scale. We change the title to: Basal freeze-on generates complex ice-sheet stratigraphy.
We adjusted the abstract on new L11-14 too, so that it reads: Here, we simulate their genesis bythe process of basal freeze-on using numerical ice-flow modelling and analyseassess the transient evolution of the emerging ice-plume and the surrounding ice-layer structure as a function of both freeze-on rate and ice flux and freeze-on rates.
Reviewer #3 ( Suppl. Table 1. It is not obvious to me that there is a spatial correlation. From Fig. 2, areas of elevated freeze-on index occur in numerous areas. In some areas, they coincide with subglacial water paths and plumes, and in other areas they do not (see also attached figure). The attempt to correlate plume locations and freeze-on index in Suppl. Table 1 does not shed light on this. Considering that the freeze-on index is high in so many parts of the ice sheet, I would expect most streamlines to cross such an area at some point.
To clarify the issue of spatial correlation; one needs to be aware that Φ only represents the geometrical requirement for freeze-on based on the pressure-induced depression of the melting-point (PiDMP) such as the relationship between the adverse bed, where water can potentially flow uphill, and the surface slope. Hence, it is a purely geometrical quantity. A steeper surface slope allows for a steeper adverse bed slope, which results in a larger freeze-on index Φ. Therefore, the index depends on the bed and surface topography, and it is a property of this index to be found favourable to PiDMP in many places.
A further point is that in order to have freeze-on water is needed. For a given water flux Q w we can calculate the amount of freeze-onḟ for each location by multiplying Φ, which is non-dimensional, with the water flux, the ice density, gravitational acceleration, and finally dividing it by the volumetric latent heat of freezing.
As we do not have direct observations of where and how much water is flowing under the ice sheet, we use Φ in our study to visualise the locations where freeze-on could possibly happen under the right conditions. It is not our aim to use Φ as a predictor for freeze-on plumes.
Next we consider the availability and quantity of water flux, which means we also need to consider the ice flux. In regions of low ice flux, freeze-on rates need not to be big to produce a visible plume while in regions of large ice flux the freeze-on rate needs to be considerable greater to obtain the same plume height.
As water paths have only been modelled and not observed so far, we assume, as an approximation, that water is everywhere and calculate the path it would take. Hence, using this method we are able to visualise the major water transport paths -the regions where we would expect potentially available Furthermore we need to be aware that the accuracy of Φ is related to the uncertainties in the surface and bed dataset. This refers in particular to the bed topography, which shows the largest spatial error as it is interpolated from RES data, where in regions of sparse radar data the vertical error in elevation/thickness exceeds 600 m (see Morlighem et al. (2017) (Supplementary Fig. 3) 14,15 ). In regions of sparse data the error along a single radar transect is between 120 − 200 m and increases away from the transect to 300 − 400 m over a distance between 5 − 10 km. This error does not represent the error in bed slope, which varies much more locally, but gives a general indication where the calculation of Φ is expected to be less accurate. Note that Φ and the water paths are consistent since the same dataset for bed and surface is used, so that with each improved topography dataset the estimates of Φ and water paths become more reliable.
Accounting for the overall error and considering the main water paths, the likely frozen areas (as well as ice flux when looking at the plume distribution) a distinct spatial correlation is visible. Some plumes are visible in areas where the likelihood map suggest a frozen bed but where various studies suggest a thawed bed. Plumes are not observed in regions where there is total agreement for a frozen bed between the studies (e.g. area G3 Response Fig. 1 rates, respectively (see Figure 5a 13 ). Compared with the eight 3-D thermomechanical models (see Figure   3 13 ) almost all models (> 90%) predict a frozen bed at the location of the green circles (except G1 and G5). Whereas for the red circles less, than half (42%) of the 3-D thermomechanical models predict a thawed bed while other sources (e.g. NorthGRIP borehole 6 , radar 4 ) observed basal melt in this region.
After we removed the threshold of very low slope from the calculation of Φ, the interior shows many more regions of Φ also within the red circles. This change in coverage leads to a higher match between observed plumes and freeze-on area (see amended Table 1). For the new calculation of the random fit, we selected the points randomly within the area covering the observed plumes and increased the number of points to 10, 000.
We undertook the following changes to improve the visualisation of the correlation: (i) Figure 2 has been amended by: using the re-calculated Φ with no artificial threshold as explained above.
choosing a coarser seed grid (45 km instead of 12 km) using a 5 km step size to visualise the main water path more clearly.
using the complete data set for plumes (41 small plumes added) as stated above.
showing the surface-velocity vectors on a 150 km coarse grid.  Furthermore, we clarify in the results section that basal freeze-on is a geometrical quantity and and make our statement in the sentence beginning on L51 less assertive owing to the uncertainties in the available datasets, and changed the sentence to highlight that the general plume distribution matches the major water paths. The paragraph now reads new L66-78 (new L78-83 shown above for Reviewer #1) now reads: Although in Greenland the relief of the subglacial topography is less pronounced than in Antarctica, we  (20)) that accounts for the effect of the bed-incline (relative to surface slope) on freeze-on triggered by PiDMP. Analysing the spatial distribution of internal ice-layer structures with regard to Φ and water flow-path, while considering the likelihood of a frozen bed, suggestsreveals that in generallarge plumes are located along the main subglacial water paths crossing areas of elevated Φ (Fig. 2, and Supplementary Figs. 2, 3 and Tables 1, 2; see Methods section 5.1 and 5.5; note, uncertainties are largest in regions of sparse radar data 15 which affects both Φ and the computed water path). According to theory, PiDMP depends dominantly on basal water-flux and topography, whereas the thermal state of the ice is of secondary importance. Thus, Oour spatial analysis thus supports the hypothesis of basal freeze-on due to ascending water for Greenland.
We expanded both methods section 4.2 (new L504-542) and 5.5 (new L597-618) on the uncertainties and clarified the geometrical aspect of the freeze-on index Φ. They now read:

Basal freeze-on index.
In order to spatially compare on an ice-sheet wide scale the potential freezeon purely related to the PiDMP mechanism, we define the basal freeze-on index Φ as an entirely geometrical quantity, which reflects the basal topography relative to the surface topography. Since the mechanism of PiDMP requires water flux along an adverse bed slope, we only take the heat terms G p and G w in equation (15) that depend on the water flux Q w and the bed-to-surface-slope ratio S to calculate Φ. The heat from the sum of the neglected heat sources G g + G s is comparable to the difference of G p − G w calculated for small water fluxes (i.e. Q w = 0.01 is m 3 s −1 per metre width) and ratios S >= −5, even for increased ice flow (i.e. 100 m a −1 ). For large water fluxes (i.e. Q w = 0.11 m 3 s −1 per metre width) however, the difference G p − G w is one to two orders of magnitude larger (Supplementary Fig. 6a). By neglecting the two heat sources G g and G s , we implicitly assume that they are being balanced by heat conduction towards the cooler ice-sheet surface as well as advection of colder surface-ice along flow.
To obtain the freeze-on index Φ as a non-dimensional quantitiy we rewrite equation (15) using the terms G p and G w only (ignoring G g and G s ) and introducing equations (18) and (17), respectively, to obtaiṅ and define Φ as where Φ now depends largely on bed and surface slopetopography and scales with Q w ( Supplementary   Figs. 6). The quantity G w is the heat generated by sub-glacial water flow, while G p is the heat needed to maintain the water at pressure-melt point. For a given water flux Q w the freeze-on rate is the result of a simple multiplication between Q w , the ice density ρ i , the gravitational acceleration g, the inverse of the volumetric latent heat of freezing 1/L, and Φḟ Comparing the non-dimensional heat termsG w andG p (Supplementary Fig. 6b; equations (17) and (18), respectively) for a range of −2 ≥ S ≥ −11 shows that the heat taken up byG p to warm the water increases with decreasing S, while the heat emitted byG w , due to water flow, decreases towards zero.
As the increase inG p is greaterstronger than the decrease inG w for decreasing S, the heat differencẽ G p −G w is increasing, leading to increased freeze-on in orderso as to balance the heat termsheats ( Supplementary Fig. 6bb).
The calculation of the freeze-on index Φ depends highly on the bed and surface slope α b and α s . However, the interpolation of the RES data to obtain the bed topography leads to large vertical errors in bed elevation exceeding 600 m in regions of sparse radar data (see Supplementary Fig. 3 in Morlighem et al. (2017) 15 ). A minimum error of 200 m allows us to nearly continously visualise all used radar transects leading to a minimum transect width of ≈ 4 km in regions of sparse data. Increasing the minimum transect width to 5 and 10 km leads to a maximum error of 225 and 300 m, respectively. Note that Φ and the water path are consistent as the same dataset for bed and surface is used, so that with each improved topography dataset the Φ and water paths become more reliable.

Estimating correlation between freeze-on index and plumes.
To quantify the relation between plumes and areas with freeze-on index areas, we follow streamlines along the reversed (surface) iceflow gradient starting at each plume, and determineevaluate if they cross areas of freeze-on index Φ.
Further, we compare the results between the mapped plume sets (≥ H/3 1/3H and < H/3 1/3H) with the averaged results for 10, 000 randomly chosen locations. For the reversed streamlines calculated on a 1050 m grid, as used in this paper, we find the best match between plumes and Φ areas for the large plumes, followed by the small plumes (result reduced by ≈ 10%) and the random set, which are reduced by up to 8% and a further 9% a result by ≈ 30 − 20% lower, respectively (see Supplementary   Tables 1 and 2). By comparing this difference between the observed and random data set relative to the total random mismatch (the percentage of plumes that are predicted to have no freeze-on area), we find that the observed data results in a distinctively better match than the random data set ( Table 2; e.g. the mismatch is reduced by > 35% and > 17% for plumes with heights ≥ H/3 and < H/3, respectively). Consequently, we argue that the plume distribution is related to the pattern in Φ, since if in the observational data set the plumes are randomly distributed we would expect a similar mismatch between the data sets.
Calculating reverse streamlines on a 150 m grid 16  ing plumes (Fig. 5). An exception occurs for very large water fluxes and along-flow series of freeze-on areas, which leads to stacking of the advected plumes (Fig. 5b)." Observations from Petermann basin (e.g.  indicate that the plumes are increasing in size downstream. The relative height might behave similarly though in absolute height the plumes in the fast flow regions are smaller than in the interior -also visible in Figure 5a,b compared to 5d; also compare Suppl. Fig. 1 with Response Figure 2. With ice flux in a converging zone, the ice flux becomes larger with both meteoric ice flux M and ice flux F originating from freeze-on contributing towards the total ice flux. For the case of uniform freezeon over the convergence area the plume height is expected to be the same within the convergence as the ice-flux ratio R hasn't changed. However, for a reduction/increase in ice-flux ratio R the converged plumes are expected to fall/rise, respectively.
In summary, our modelling and suggested mechanism of freeze-on is fully consistent with the Petermann Response Figure   In Fig. 2, there is a large area with elevated freeze-on index but the plumes are downstream of this area.
This seems to be in direct contrast with the modelling results (e.g. Fig. 5b where the plumes decrease in size downstream). Are the plumes observed in Petermann a result of "very large water fluxes"? And in this case, how large? Can the model reproduce the observations from Petermann?
Looking at the area with elevated freeze-on index we realised that we 'lost' some of the small plumes originally displayed in the Supplementary Fig. 1 of the originally submitted paper. This has now been amended so that all plumes from the original dataset ( Fig. 1 and Supplementary Fig. 1) and the dataset added during revision 1 are used in the plots.
The relatively wide strip with elevated Φ upstream of Petermann glacier ( Figure 2 in main text) is again in an area of a likely frozen bed. The distribution of the water flow paths suggests that only a few main paths (2)(3) are crossing this area, coinciding with the occurrence of plumes. This fits nicely with our suggested mechanism of basal freeze-on.
With the previously missing 'small' plumes added it is evident that the plumes start where the freeze-on index is strong, with plumes first being small and increasing further downstream. As the transects are across flow the longitudinal extent of the plume is not visible, but one can see that downstream the observed plumes on the across transects become bigger. This is consistent with the modelling results in Lines 199-202: "Away from the divide, the first large plumes (> 1/3H) start to form in regions where the basal water flux Qw is increasing, with the total ice flux Q and freeze-on index Φ still being generally small. Further along, towards the margin both Φ and Qw tend to increase (while Q is still moderate) and allow for larger plume heights." Based on the suppl. Fig. 2, I would argue that it is just as likely that the plumes form in areas where there is a transition from frozen to thawed conditions.
It is possible that transitional areas between frozen and thawed conditions are a good environment for plume growth (especially < H/3, as in such locations a limited amount of water would be expected).
However, for plume heights > H/3 this argument is not as strong as our observations made in context with Φ, water flux and ice flux, as in this case we would expect large plumes everywhere along this transition. By plotting the data with the main water paths it is generally visible that large plumes are located at the start or along these water systems (new Fig. 2 ).
We changed the sentence by adding that the regions likely correspond to basal melt water converging, so that it now reads (new L228-230): Away from the divide, the first large plumes (> H/3 1/3H) start to form in regions where the basal melt water is likely convergent, leading to an increasingthe basal water flux Q w is increasing, with the total ice flux Q and freeze-on index Φ still being generally small.

#2 Amount of basal meltwater and freeze-on rates
Line 154: Reference no. 27 and 28.
Both references here are modelling studies and do not provide direct evidence of subglacial water or major pathways.
We agree entirely with the Reviewer. The text does highlight this, as we wrote the following in the sentence ending on line 154 (now new L186-189): 'However, recent modelling studies suggest the existence of major meltwater pathways, subglacial lakes and, related to the mechanism described as the pressure-induced depression of the melting-point (PiDMP), widely occurring basal freeze-on 7,12 . ' Line 167-168: "It is readily conceivable that locally conditions arise, which lead to freeze-on rates comparable to the ones used in our experiments." Is it? This statement comes across as unsubstantiated. Is the slope ratio reasonable for the areas where the plumes are observed? Are the values for the water flux reasonable? This needs to be stated explicitly.
Also, there is a difference between conditions such as these arising locally, and to them being present on a large scale across the entire ice sheet as evidenced by the extensive presence of the plumes.
Although Φ is observed on a large scale over the ice sheet it is a local phenomena that depends (as explained earlier) on the bed and surface slopes locally (grid size). The smaller the grid used to calculate the slopes, the larger the range in local slopes (flatter as well as steeper).
We changed the sentence to refer to the grid size (new L201-203): As the calculated bed and surface slope depend heavily on the used grid size, itIt is readily conceivable that locally conditions arise in Φ and Q w that lead to freeze-on rates comparable to the ones used in our experiments (Fig. 7).
Lines 251-252: "and broadly coincide with modelled freeze-on regions [28]" This statement is a bit misleading since the freeze-on rates in the study cited (Dow et al.) are at least an order of magnitude smaller than what is invoked here.
In our text we state that the regions broadly coincide and not the freeze-on rates. Dow et al. (2018) 7 state in their discussion that they find subglacial water accumulating in troughs feeding NEGIS and Petermann Glacier (as well as Jakobshavn Isbrae), which allows for more supercooling freeze-on to occur (however with low rates).
The freeze-on rates are expected to differ with our study as the grid chosen by Dow Supplementary Fig. 1)). We used a constant water flux of Q w = 0.11 to multiply with Φ in order to get the freeze-on ratesḟ over these two regions for four different grid sizes (150, 300, 1, 050 and 4, 950 m). The maximum value in freeze-on rate for the entire domain (freeze-on rate calculated for the given Q w for the maximum Φ found in the domain for each grid size) is 6 times larger for the smallest grid than for the largest grid (see Table 2. However the location is likely not the same. Comparing at two locations we see from the Table 2 that the values vary from site to site for different grid sizes. From this example one can see that the freeze-on rates depend locally very much on the grid size. On a coarser grid, the slopes get averaged out, meaning the big picture is the same but locally the values differ, so that the freeze-on rates calculated for a 5 or 1 km grid are substantially different. Also, the water flux is expected to change locally to with a different grid size -as water paths would change -and affect the freeze-on rates (see Figure 7 in main text).
Moreover, the study by   7 refers in general to averaged freeze-on rates, mentioning Qw, the ice density ρ i , gravitational acceleration g, and dividing it by the volumetric latent heat of freezing L for the whole modelling domain using four different grid sizes. The maximum value found for the whole domain as well as the maximum value found in two specific regions with a patch of Φ are shown (upstream Petermann (centered at −80 km easting and −1240 km northing) and the region with the highlighted plume from the paper (Supplementary Fig. 1). We added a sentence to the discussion (new L298-300); Discrepancies in freeze-on rates between models and reality can be explained by differences in topographic resolution, affecting both the water pathways and the range in bed and surface slopes, with a coarser grid levelling out the extremes in slope.
Lines 262-263: "Therefore, statements about the mismatch between plume size and freeze-on rates [28] are not valid without considering as well the effect on the relationship between plume height and freezeon rate both from the flux relationship and rheology contrast." I agree that looking at the freeze-on rate in itself is not adequate when examining plume heights. However, from the figures it seems that 10*softer ice roughly leads to a doubling in plume height. Since the freeze-on rates are one to two orders of magnitude larger than the rates from Dow et al., using the freeze-on rates from the latter study would never lead to a plume of that height, regardless of whether We added a section in the Methods (new section 3.1) and referred to it in the main text (L141-149 and L196-197)) so that the statement on Lines 262-263 make more sense to the reader.
The paragraph new L175-183 is slightly changed to: The flux-dependence illustrated has crucial implications for areas of low surface-accumulation rates (e.g. East Antarctica or Northeast Greenland), since where meteoric ice flux M is small, even a low freeze-on flux F is able to produce relatively large plumes (see Methods section 3.1). This effect is amplified in regions with small surface accumulation areas, which further reduceing the meteoric ice flux M ; such regions, which exist for instance near ice divides or wherediverging flow is divergentregions, typically upstream ofabove rising beds 8,17 . Downstream, plume heights generally decrease owingdue to the increasing meteoric ice flux M , which is a consequence of both contributing accumulation area and flow-convergence. In addition, basal sliding becomes more important downstream, affecting the distribution of horizontal velocity with depth, and thus further decreasesing the plume height (Fig. 6).
And Lines 196-197 are changed to new L232-233: As the ice flux is low, only low water fluxes are needed to produce the required freeze-on rates (a few millimeters to centimeters per year; see Methods section 3.1).

#3 Frozen basal conditions
Lines 194-197: "Near the ice divide the observed plumes are usually smaller than 1/10H, in the range of a few 100 meters. For these relative small plume heights the flux ratio R too is small when assuming flow purely by internal deformation. As the ice flux is low, only low water fluxes are needed to produce the required freeze-on rates (a few millimeters to centimeters)." Near the ice divide, radar observations indicate that the bed is likely frozen (as shown in Suppl. Fig.   2), thus the subglacial water flux must be zero and therefore also the freeze-on rates. It is a significant weakness in this manuscript that the presence of the plumes in the interior of the ice sheet cannot be explained by the freeze-on process but this weakness is not mentioned or discussed.
Many of the changes made so far to the revised manuscript together with answers given in this letter contribute to the above comment. Nevertheless, we reply to this comment with some more detail: As mentioned above we improved Figure 2 by adding the likely thermal state for a frozen bed (together with the outline for low ice flux; see answer to comment #1 by Reviewer #3) and moved the original figure 2 to the Supplement (now Supplementary Fig. 2).
However, below follow some more arguments why we cannot exclude the presence of water and therefore plumes near the ice divide.  13 highlight that their synthesis of thermomechanical models and remote inferences agree in general for the Northeast Greenland Ice stream (thawed) and the west facing slopes of the central ice divides (frozen), but that e.g. the borehole-observed basal thermal state near NGRIP (thawed bed) is rarely represented by either the models or remote inferences and leaving approx. one third of the ice sheet in need of more observations. Therefore, we cannot exclude that some water paths in regions of a likely frozen bed exist, but rather assume that the likelihood of water is smaller. As an example Rogozhina et al. (2016) 19 predict a high geothermal heat flux in the region south of NEEM and to the southeast, so that in these regions where a likely frozen bed is predicted 13 a thawed bed instead might be likely (Fig. 3a in Rogozhina et al. (2016) 19 . This would also fit with our observations of existing plumes. We could argue that if these observed near basal structures are due to freeze-on we can use the plumes to inform on possible water paths beneath the ice sheet and in consequences make inferences about the geothermal flux in particular.
Looking at the two red circles along the divide, highlighted by the Reviewer #3 in the review attachment, we find that they are in regions with very low ice flux and very flat surfaces. The likelihood of thermal state predicts this region to be a thermal boundary, but from the borehole at NGRIP and from radio echo sounding 4,13 it is known that in this region there is melting at the ice-sheet bed. Furthermore, the divide region has a very flat surface which easily allows for stronger negative bed-to-surface-slope ratio (S), with freeze-on rates approaching the centimetre per year range even for small water fluxes (see Figure   7a in the main text). Owing to the very low ice flux, these freeze-on rates are enough to lead to observed plume heights between 100 to 300 m. The red circled area further to the east is still in an area with low ice flux (balance velocity ≈ 10 m a −1 , however with a steeper surface slope (see Figure 1 in main text) and likely more water flowing at the base, both allowing for larger freeze-on rates (larger Φ) compensating for the increased ice flux.
We added a paragraph to the introduction to highlight how poorly constrained the basal temperature of the Greenland ice sheet is (new L57-63).
The basal temperature of the Greenland ice sheet is only poorly constrained, as we do not know the geothermal heat flux or how it varies on a kilometre to 100 km scale. Analysis combining independent data are needed to obtain reconstructions for the basal temperature 13,19 . Neither do we know the basal topography well enough to evaluate if there are warm-based deeper parts surrounded by frozen areas 5,15 .
In addition does the community's limited understanding of basal hydrology prevent us from predicting its thermodynamics consequences, in particular answering whether it is possible for melt-water channels to advance through freeze-on areas so as to connect isolated melt patches.

20
The latest revisions contribute further to over the clarity of the results and manuscript. I support its publication.
In my opinion, the additional comments raised by Reviewer 3's have been adequately addressed.
In particular: -The author can clearly justify how their work differs and improves upon other related studies.
-The amendments brought to the new Figure 2, including the additional plumes and removal of artificial threshold allows the authors to further support the existence of a spatial correlation between plumes, freeze-on index and water paths. It is true that uncertainties exist, and the authors now acknowledge this more explicitly in the text.
-Response Figure 2, supported by comments, reconciles the author's work with observations from Petermann basin -There are further details on the distribution and amount of freezing, involving mainly their sensitivity to grid resolution.
-The last point concerns the possibility of freeze-on / basal melting in the interior. The main argument brought by the author relates to the overall poor knowledge of the basal conditions in Greenland. They added a new paragraph now discussing this point openly, and as such, I believe have also adequately addressed the last concern.