Summary:

This paper is a modeling study of subglacial sediment transport. Subglacial sediment production and transport is important as a major landscape construction process in glacial-deglacial settings, and can greatly impact downstream rivers’ water quality. Whereas  physics-based glacial ice flow models are now sophisticated enough to be merged into global climate models, understanding of the role of basal erosion and subglacial sediment transport interactions is much less understood. The authors are pioneering a reduced-complexity numerical model of glacial processes and subglacial hydrology.

The study implements a two-dimensional subglacial sediment transport model using previously formulated subglacial sediment transport and bedrock erosional processes presented in Delaney et al 2019. A new, be it simple, planview flow and sediment routing scheme is added to the model to transports sediment down-glacier based upon the hydraulic potential gradient. It is in this coupled component that the novelty of this paper lies!

However, findings from some of the synthetic model case studies (most notably the ice sheet case) are hard to interpret/generalize, as process outcomes depend on the simplifications and boundary conditions. I recommend improvements to the study design and results sections. These changes involve omission and clarification, they do not change the model design perse or require entirely new simulations, so I would call them moderate revisions. I detail suggestions in the notes below.

Major comments

The study implements a two-dimensional subglacial sediment transport model using previously formulated subglacial sediment transport and bedrock erosional processes presented in Delaney et al 2019. A new, be it simple, planview flow and sediment routing scheme is added to the model to transports sediment down-glacier based upon the hydraulic potential gradient. It is in this component that the novelty of this paper lies.

Three cases are used to demonstrate model input and behavior. These synthetic cases are inspired by the Subglacial hydrology Model Intercomparison Project (SMIP), which is a sound applaudable approach! The aim is to use the synthetic test cases to show the model’s ability to reproduce known processes and find new insights into the spatially distributed processes responsible for subglacial sediment dynamics. However, I think that there may be improvements to the study design and results sections needed.

For the synthetic ice sheet case, this study finds that sediment discharge in the model run decreases after the climate warms and reaches a stable regime. The decrease occurs due to sediment exhaustion from increased water discharge, and associated sediment transport, which removed till that was unable to be transported in a cooler climate with less available meltwater able to transport sediment.

I find this model behavior problematic to generalize: it would be an interesting finding, depletion of sediment flux, yet the outcome is entirely controlled by the model assumptions on initial till height and bedrock sediment production. So, it is hard to draw any conclusions on this process being of relevance?

I think the paper would be strengthened by omitting this case study, it oversimplifies the ice sheet system by too much.

Results on the effect of spatially-distributed fluxes seem of more of importance, and this effect is more pronounced in the alpine glacier case. The alpine glacier case recreates high sediment concentrations in early season, and his is explained by spatial variability in the till distribution and sediment transport access to these sediment patches. This effect has been observed in proglacial streams.

It seems that there is a large process parametrization discrepancy between the general bedrock sliding law applied (with no grainsize dependence and validated over glacial-deglacial scales) and the heavily grainsize dependent sediment transport law (Engelund-Hansen). A note of caution should probably already be added into the results section, and refer to the later discussion section on this topic.

On the structure of the paper:

Title: I recommend simplifying the title.

Suggestions: “Bedrock erosion and sediment transport variations across a glacier bed controlled by glacier behavior and hydrology” or “Modeling of the spatially distributed nature of subglacial sediment erosion and transport dynamics”?

Introduction: can be tightened. Please read carefully through and omit some of the sections that are repetitive or wander.

Model description: I do appreciate the review of the hydraulic model, even if it previously been well described in Delaney, 2019, but is needed here to have this paper be a stand-alone contribution.

I found the order in Section 2.2 non-intuitive. Indeed, it is important to detail the Exner equation approach first, but  then begin with bedrock erosion parametrization, as that is the term that produces sediment, and then describe sediment transport.

One suggestion to make each of the source terms and in-out fluxes clear is to add a diagram of a mass fluxes in the Exner equation, and label the processes. Like Fig 1 in Paola and Voller, but then made specific for this special implementation.

Line-by-line minor comments:

line 6: the concept of 'sediment connectivity' needs explanation or definition before being invoked as an important control. Either explain within the abstract, or omit.

Line 14: this sentence is unclear."We find that sediment grainsize plays an important role. Smaller sediment sizes....

Line 24: possibly add this paper by Dongfeng Li 2021, although the increasing sediment loads are not just attributed to glacier melt but also due to permafrost thaw and rainfall change. Li, D., Lu, X, Overeem, I., Walling, D.,  Syvitski, J.,  Kettner, A.J.,  Bookhagen B.,  Zhou, Y., Zhang, T.,  2021). Exceptional increases in fluvial sediment fluxes in a warmer and wetter High Mountain Asia. Science, 10.1126/science.abi9649

Line 30: replace with: 'are limiting nutrients in the oceanic ecosystem’

Line 45: add flow after …water

Line 51: although if there is little sediment embedded there is little abrasion!

Line 90: please clarify this approach: is there a single hydraulic diameter and single associated water discharge for the entire distributed drainage system? Or dothese properties Q* and Dh get assigned/calculated for individual drainage channels?

Line 94: this selection of representative discharge seems really difficult and perhaps arbitrary? how do you decide ahead of time what quantile (or did you mean quartile?).

Line 103: please spell out R-channel at its first occurrence.

Line 124: I am not sure about the limit, H-lim, at 10 cm, this seems really arbitrary. Sediment transport in a pressurized pipe flow can probably easily scour bedforms to 4-5 times that depth almost instantaneously? I understand that perhaps there is little data to constrain this parameter, but it may be prudent to explore whether the model is sensitive to this setting.

line 165: this switch in describing code is a bit out of style, perhaps better to describe this in the code documentation as opposed to command lines in the paper.

line 177 Reference is oddly formatted, please make conform journal requirements

line 184 First explain why this would be a non-stiff problem, and then state that this solver is apprporiate for non-stiff problems.

Line 202-208 this section has several ‘disclaimers’ on your assumptions that would be better suited for the discussion of model limitations later in the paper.

Line 245: I am not intimately familiar with SHMIP, but does the rate of temperature offset originate from those setups? Do I understand it right that the total warming scenario is an added 15 degrees C to diurnal amplitude? Over 30 years? That warming rate seems really abrupt, and unprecedented even under the most catastrophic warming scenarios?

Line 255: line 255 replace till height by thickness?

Line 372. Repeat for careless readers, what are the last three time spans?

Line 374. Isn’t one of the important changes of this approach that you do have a way to access different patches of the glacier bed?

Could it be model underestimation of till thickness and erosion instead?

Figures

Figure 2: This figure is really helpful in providing a feel for the dimensions of the subglacial drainage network , flow velocity magnitude (small), and spatial distribution of discharge.

Figure 3. Recommend this figure to go to an appendix. It is not too well explained in the text and seems important for a manual of the model, not for the scientific findings.

Figure 5: Perhaps improve this figure by changing the aspect ratio of these figures, at present they are hard to read.

May help to add 2 initial panels with icesheet topography map and ice flow velocity? And then have the 'maps' of the basal conditions.

Figure 7

May help to add 2 initial panels with icesheet topography map and ice flow velocity? And then have the 'maps' of the basal conditions.

Perhaps improve this figure by changing the aspect ratio of these figures, at present they are hard to read.

Great to see this as animations.

Figure 9. Remove ‘likely ‘ from the caption.

Figure 10.Add to the figure caption what parameters are modeled vs observed.

As a preliminary note, I would like to state that F. Herman, who is last  author of this manuscript, was also last author of a paper where I was involved (Prasicek et al., EPSL, 2020). However, I never met the first and second author of this manuscript in person, and there is no further connection. So there is no conflict of interest in any direction. As a second note, I would like to mention that I have not worked in the field of glacial erosion and sediment transport for a long time. So I may be not familiar with some of the concepts used in modeling glacial landform evolution. In turn, I have worked with fluvial models and hydrogeological models for a long time and believe that my mathematical background might be helpful here.



Owing to my background, my review focuses on the methods part, while I wrote only a few remarks on the results and discussion sections.



The topic of subglacial sediment transport is important in the field of glacial landform evolution, and it is perhaps the component with the biggest gaps in knowledge. In my own work, I used simple analogies from fluvial sediment transport with almost no validation by real-world data or by more elaborate models. So I find this topic very interesting, but I am quite sure that it is interesting and potentially important for the glacial erosion community.



On the other hand, the description of the theory and the numerical part falls much behind my expectation on the group of authors and also much behind the 2019 JGR paper about the 1D version. At some places, it even looks as if it was wrong. So think think that considerable parts of the theory need to be rewritten.

(1) Sect. 2.2

The main part where I am not immediately convinced that is is correct is Sect. 2.2 about sediment transport, starting from the balance (Exner) equation (Eq. 4). In its genuine form, such an equation would be written in 2D in terms of a sediment flux per unit width (m^2/s) instead of the sediment flux (m^3/s). In the 1D version (Delaney et al., JGR, 2019, Eq. 9), the sediment flux is used in combination with a channel width w. However, the version introduced here uses a length scale l that describes a "characteristic length-scale for sediment mobilization, over which sediment mobilization adjusts to sediment transport conditions." This length scale is a longitudinal length scale (in flow direction), while the formulation of the balance equation in terms of the flux requires a length scale perpendicular to the flow direction (such as the channel width). So a am quite sure that Eq. 4 is not correct in the form it is written, but I cannot assess whether it is correct in the implementation and whether it affects the results in case it is not correct.



There also seems to be a problem with the physical dimensions in Eq. 5 (beyond that the divergence of the sediment flux at the left-hand size should not be called "sediment discharge"). If dot(m)_t is indeed an erosion rate (m/s) as defined earlier, it is not consistent with the other properties, which are fluxes per length (m^2/s).

(2) Sect. 2.3

While I am confident that the numerical implementation is sound, I am not happy about the way it is described in this section. To be honest, I found it even more confusing than enlightening. The cited work by Bovy et al. (2016) used a standard finite-volume discretization of the fluxes, while the algorithm proposed by Braun and Willett (2013) used a single-flow-direction (D8) scheme (if I am not wrong). I read that you use a multi-flow-direction scheme, but it is not clear what the difference toward typical continuum schemes (finite volume) is or whether it is the same. Instead of (or in addition to) mentioning functions and software packages, I would ask you to state the equations that are finally solved.

(a) How exactly is the flux (water/sediment) distributed among the available directions? It looks as if this changes trough time, but I did not get it completely.

(b) What does "integrating Equation 1" (line 166) exactly mean here? Typically, "integrating" a differential equation is somehow one-dimensional or upstream in a tree, like the algorithm made popular by Braun and Willett (2013). However, if there are multiple flow directions, upstream paths may meet again, so that you would end up at different values of phi by integrating over different paths. I see that you are solving a system of equations, but I cannot see how.

(c) It would be good to state the balance equation for the water and for the sediment for each grid cell explicitly and in such a way that it becomes clear which system of equations finally has to be solved.

(d) Also about the discretization: I guess you are assuming one single conduit of a given hydraulic radius for each grid cell. This may be questionable, but of course not necessarily bad. A similar assumption is typically made in karst evolution models. It should be discussed at least briefly. And is this conduit directed, or will it change its direction if phi changes?



There are also some more parts in the model description that should be clarified, where I think this can be done quite easily.

(3) Lines 104-105: Q*_w or Q_w? And I did not understand how the water pressures can increase to unreasonable values if the flux rapidly increases. I can imagine that this happens if D_h is small for a single-flow-direction model, but why do the alternative flow directions not help here? (4) Line 137: Here you use "capital S" (the cross section area), which is definitely different flow "lowercase s" used in Eq. 1 (probably a nondimensional factor).

(4) Line 137: Here you use "capital S" (the cross section area), which is definitely different flow "lowercase s" used in Eq. 1 (probably a nondimensional factor).

(5) Line 144: As far as I know, the lower bound l_er = 2/3 mentioned here was obtained from a worldwide comparison of glaciers under different conditions. So I am not sure whether this low value is relevant for Eq. 10.

(6) Lines 144-149: Some other studies use the relation from Eq. 11 for the deformation velocity and a similar relation with an exponent n-1 instead of n+1 for the sliding velocity u_b. This relation predicts that deformation becomes more and more relevant if the thickness of the ice layer z_s-z_b increases. Using this relation in my own work, I was even told by reviewers that the sliding velocity cannot be predicted from the deformation velocity at all. I do not share this point of view, but I think it should be discussed why even a weaker assumption than the relation typically used is employed here.



I am not familiar to the test cases used in Sect. 3, so that I just list some points that came into my mind when reading it.

(7) Line 239: Why is there a need to increase DeltaT, and what is the value of DeltaT used here?

(8) Line 246: I did not get the point why H_max grows.

(9) Line 285: What is the meaning of the -5 in Eq. 16, given that DeltaT is an (adjustable) temperature offset?

(10) Lines 339-342: Nice numbers, but you somehow let the chance to analyze your model more thoroughly pass by. From Eqs. 7 and 8, we see that the transport capacity Q_sc is proportional to v^5 and thus also to Q_w^5. this is what you put into your model. So you could integrate Q_w^5 over, say 1 year periods, and look how well your sediment output correlates to this integral. If it correlated perfectly, then your sediment output was just what you put into your model equations, and everything else would be unimportant. I think it will not correlate perfectly, so that you can discuss which of the components of your model is important. However, this is just an idea how you could sharpen the discussion.

Please do not get me wrong -- it is definitely not my intention to tear down your work. However, there are already many modeling papers in the literature where readers cannot see clearly enough what was done and reproduce the basic ideas. As modelers, it is our task to act against this tendency. I am quite sure that this piece of work will finally end up at a very good paper, but you will have to spent some more work on it.



Best regards,

Stefan Hergarten

Report attached