The spatially distributed nature of subglacial sediment dynamics: using a numerical model to quantify sediment transport and bedrock erosion across a glacier bed in response to glacier behavior and hydrology

Abstract. In addition to ice and water, glaciers expel sediment. As a result, changing glacier dynamics and melt will result in changes to glacier erosion and sediment discharge, which can impact the landscape surrounding retreating glaciers, as well as communities and ecosystems downstream. To date, the available models of subglacial sediment transport on the sub-hourly to decadal-scale exist in one dimension, usually along a glacier's flow line. Such models have proven useful in describing the formation of landforms, the impact of sediment transport on glacier dynamics, the interactions between climate, glacier dynamics, and erosion. However, because of the large role of sediment connectivity in determining sediment discharge, the geoscience community needs modeling frameworks that describe subglacial sediment discharge in two spatial dimensions over time. Here, we present SUGSET_2D, a numerical model that evolves a two-dimensional subglacial till layer in response to the erosion of bedrock and changing sediment transport conditions below the glacier. Experiments employed on test cases of synthetic ice sheets and alpine glaciers demonstrate the heterogeneity in sediment transport across a glacier's bed. Furthermore, the experiments show the non-linear increase in sediment discharge following increased glacier melt. Lastly, we apply the model to Griesgletscher in the Swiss Alps where we use a parameter search to test model outputs against annual observations of sediment discharge measured from the glacier. The model captures the glacier's inter-annual variability and quantities of sediment discharge. Furthermore, the model's capacity to represent the data depends greatly on the grain size of sediment. Smaller sediment sizes allow sediment transport to occur in regions of the bed with reduced water flow and channel size, effectively increasing sediment connectivity into the main channels. Model outputs from the three test-cases together show the importance of considering heterogeneities in water discharge and sediment availability in two dimensions.


I'm very surprised to see the exponent 5/2 in equation (7). Most of these bedload transport laws have 3/2. I don't have the Engelund-Hansen report, but in his 1970 JFM paper, he uses Meyer-Peter/Müller (and doesn't reference this report). Ah, I see reading on (140) that there is a reason for this, as it supposedly includes suspended sediment. Of course, a proper treatment of suspended sediment then requires an evolution equation for the suspended sediment concentration. It seems to me that if you go to the trouble to include bed erosion to the bedload layer, then it is logically commensurate to include suspended load concentration, at least in some fashion. But again, we must be thinking of streams, and then such streams do not cover the glacier bed; how should sediment transport to the streams be modelled? The statement "the continuous nature of the relationship improves the model stability" is poor. First you pose the model, then you deal with trying to solve it. You don't decide what is in the model on the basis of what you can solve (or you shouldn't).
At equation (11), I begin to wonder what is the point of this exercise. Yes, what happens at the bed is complicated. But to choose the sliding velocity to be a fixed fraction of the shearing velocity is simply making things up. Particularly, sliding depends on the subglacial hydrology through its dependence on the effective pressure. One might argue that the emphasis here is on sediment transport, but that fundamentally depends on water flow (and also actually till deformation), and I see little point in trying to deal with one at the expense of the others, at least if the results of the model aim to be realistic.
As we come to the numerical implementation, I belatedly realise that the point of all the simplifications to the physics is that it allows a cellular model to be constructed. It reminds me a bit of the paper by Barchyn et al. (2016). I am not a big fan of this approach, which seems to me to be motivated by the wish to produce pretty pictures at the expense of doing science. A model is only as good as its formulation, and I find the modelling here to be weak in a number of points.
This point is perhaps illustrated by the comment at 371, the 'model successfully captured the inter-annual variability in sediment discharge from the Griesgletscher'. But a parameter search was used to find parameter values which worked. So can we conclude that the model is a good one? No. Does it then have predictive value? No. And, most importantly, should we expect it to be a good representation of physical process? In view of my comments above, I would have to say no.

Smaller points:
At equation 1. This looks a bit odd to me. In my way of thinking, for force balance in a channel, you would have τ l = ρgSA, where τ is the stress, l the wetted perimeter, S a suitable slope and A the cross-sectional area. So Ψ = ρgS = τ l A , and if τ = f ρu 2 and Q = Au, then Ψ = f ρQ 2 l A 3 . You can get equation (1) if l ∝ A 1/2 , as for a circular or triangular cross section; but it seems to me that equation (1) is not the basic law. For example, a wide stream has l ∝ A, approximately.
Line 105. This is unclear. D h is a length, not an area. 115. and → or, presumably.
181. The wording here suggests that a partial differential equation is being solved, but if I understand this correctly, this is not the case. (4) with (5) form a set of ordinary differential equations.