Modeling the spatially distributed nature of subglacial sediment transport and erosion

Glaciers expel sediment as they melt, in addition to ice and water. As a result, changing glacier dynamics and melt produce changes to glacier erosion and sediment discharge, which can impact the landscape surrounding retreating glaciers, as well as communities and ecosystems downstream. Currently, numerical models that transport subglacial sediment on sub-hourly to decadal scales are one-dimensional, usually along a glacier's flow line. Such models have proven useful in describing the formation of glacial landforms, the impact of sediment transport on glacier dynamics, and the interactions among climate, glacier dynamics, and erosion. However, these models omit the two-dimensional spatial distribution of sediment and its impact on sediment connectivity – the movement of sediment between its detachment in source areas and its deposition in sinks. Here, we present a numerical model that fulfills a need for predictive frameworks that describe subglacial sediment discharge in two spatial dimensions (x and y) over time. SUGSET_2D evolves a two-dimensional subglacial till layer in response to bedrock erosion and changing sediment transport conditions. Numerical experiments performed using an idealized alpine glacier illustrate the heterogeneity in sediment transport and bedrock erosion below the glacier. An increase in sediment discharge follows increased glacier melt, as has been documented in field observations and other numerical experiments. We also apply the model to a real alpine glacier, Griesgletscher in the Swiss Alps, where we compare outputs with annual measurements of sediment discharge. SUGSET_2D accurately reproduces the general quantities of sediment discharge and the year-to-year sediment discharge pattern measured at the glacier terminus. The model's ability to match the measured data depends on the tunable sediment grain size parameter, which controls subglacial sediment transport capacity. Smaller grain 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. The model provides the essential components of modeling subglacial sediment discharge on seasonal to decadal timescales and reveals the importance of including spatial heterogeneities in water discharge and sediment transport in both the x and y dimensions in evaluating sediment discharge.

1 Introduction

Increasing glacier ablation perturbs the processes through which glaciers erode bedrock and supply sediment downstream (e.g., Church and Ryder, 1972; Lane et al., 2017; Delaney and Adhikari, 2020). Changing sediment discharge from glaciers in alpine and polar landscapes impacts many social and Earth systems (Milner et al., 2017; Li et al., 2021). Thus, predictive models are needed to understand the response of these systems to glacier retreat. In alpine environments, increased sediment discharge leads to the more rapid filling of proglacial reservoirs (Thapa et al., 2005; Li et al., 2022) and abrasion of hydropower infrastructure (e.g., Felix et al., 2016). The flux of sediment from glaciers also dramatically alters alpine ecosystems (Milner et al., 2017) and occurs when high melt extends up-glacier, thus mobilizing sediment in new areas (Lane et al., 2017; Delaney and Adhikari, 2020; Li et al., 2021). In Arctic environments, increased sediment discharge can affect biogeochemical cycles given that sediments may carry phosphorus and iron (Bhatia et al., 2013; Hawkings et al., 2014). These elements are limiting nutrients in the oceanic ecosystem, so any change to sediment discharge from glaciers and ice sheets therefore alters Arctic ecosystems (Wadham et al., 2019).

Generally, two processes determine the sediment discharge below glaciers: one process adds sediment, and the other removes sediment from subglacial till layers (Fig. 1; Brinkerhoff et al., 2017; Delaney et al., 2019). Bedrock erosion adds material to the subglacial till layer. This is accomplished by quarrying, when pressure differentials on opposing sides of obstacles cause fractures to expand and rock to detach (Iverson, 1990, 2012; Alley et al., 1997; Hallet et al., 1996), and by abrasion, when debris embedded in the ice grinds bedrock as the glacier slides above (Hallet, 1979; Alley et al., 1997). Representation of these physical processes in models requires independent knowledge of a large number of parameters (see Ugelvig et al., 2018), so many researchers use empirical relationships that relate glacier sliding to glacier erosion (Humphrey and Raymond, 1994; Koppes et al., 2015; Herman et al., 2015; Cook et al., 2020). These relationships prove especially useful when applied over large temporal and spatial scales, for example, to explore the coupling of glacier erosion, climate, and tectonic uplift (e.g., Egholm et al., 2009; Prasicek et al., 2018; Herman et al., 2018; Prasicek et al., 2020; Seguinot and Delaney, 2021).

Figure 1 Cartoon of erosional and sediment transport processes considered in the model overlaid on an image of Griesgletscher in 2016. Bedrock erosion scales with sliding speed (u_s) and adds material to the till layer with thickness (H), while water (Q_w) transports sediment (Q_s) fluvially if sediment persists in that location of the glacier bed and fluvial conditions are sufficient for transport. Photo credit: Ian Delaney.

Conversely, fluvial sediment transport mobilizes material from subglacial till layers (e.g., Walder and Fowler, 1994; Ng, 2000) or it may deposit the mobilized material in the till layers under certain hydraulic conditions (e.g., Beaud et al., 2018a; Hewitt and Creyts, 2019; Delaney and Anderson, 2022). As subglacial water velocity increases, the sediment of a given grain size is transported below the glacier, and this sediment may be deposited if the water velocity slows (e.g., Meyer-Peter and Müller, 1948; Paola and Voller, 2005). Sediment mobilization ceases when no sediment is present, and the system is supply-limited (e.g., Mao et al., 2014). Thus, fluvial sediment transport depends on both the subglacial hydraulic characteristics (e.g., Walder and Fowler, 1994) and the availability of sediment at the glacier bed (e.g., Willis et al., 1996; Swift et al., 2005).

Bedrock erosion and fluvial sediment transport vary depending on the characteristics of each glacier. Bedrock erosion processes tend to dominate sediment discharge below glaciers with minimal sediment storage, large concentrations of subglacial debris entrained at the glacier bed, and steep gradients (Hallet, 1979; Humphrey and Raymond, 1994; Herman et al., 2015, 2021; Ugelvig et al., 2018). Landscape evolution models that represent glacier landscapes illustrate the dominant role of erosional processes, as opposed to sediment transport processes, over geologic timescales (Harbor et al., 1988; Herman et al., 2011; Egholm et al., 2012). Over shorter timescales of months to decades, however, fluvial sediment transport often drives sediment discharge (e.g., Delaney et al., 2018b, 2019; Perolo et al., 2018).

The development of numerical models of subglacial sediment transport to date has mainly focused on processes acting in the down-glacier (e.g., x) dimension. To date, one-dimensional models have yielded insights into the creation of eskers (Beaud et al., 2018a; Hewitt and Creyts, 2019), the formation of subglacial canals through which water flows (Walder and Fowler, 1994; Ng, 2000; Kasmalkar et al., 2019), subglacial processes in overdeepenings (Creyts et al., 2013), and the behavior of tidewater glaciers (Brinkerhoff et al., 2017). Yet spatial heterogeneities in the distribution of sediment and sediment transport capacity (largely controlled by water velocity) commonly result in the water carrying less sediment than could be transported theoretically (e.g., Delaney et al., 2018b). As a result, reducing the modeled processes to one dimension omits key processes controlling sediment dynamics because subglacial water flows through spatially distributed networks of cavities and channels across the glacier bed (e.g., Werder et al., 2013). Therefore, describing subglacial sediment transport inherently lends itself to a discretization of bedrock erosion, sediment transport, water flow, and sediment availability in both the down-glacier and transverse dimensions (e.g., x and y).

Here, we present SUGSET_2D, a two-dimensional subglacial sediment transport model that includes sediment transport and bedrock erosion processes. We implement a routing scheme that transports sediment in x and y directions based on the local hydraulic potential gradient. Synthetic cases demonstrate the model's ability to reproduce known processes and yield insights into the spatially distributed processes responsible for subglacial sediment dynamics. We also apply the model to a real alpine glacier, Griesgletscher, in Switzerland, which has a record of subglacial sediment transport from the catchment from the period 2011 to 2016. The model was run with topography data and modeled hydrology from the glacier, and measured sediment discharge data were used to validate the model. Through these experiments, we identify key processes of sediment transport from subglacial environments.

2 Model description

The model presented here implements a hydraulic model and sediment routing scheme that translates many of the underpinnings of the one-dimensional subglacial sediment transport model presented in Delaney et al. (2019) to two dimensions. We describe hydraulic and sediment transport models, explain the implemented water and sediment routing scheme, and outline its numerical implementation in two dimensions.

2.1 Hydraulic model

SUGSET_2D requires a hydraulic model as a means to route sediment and water through the subglacial environment. The hydraulic model determines the sediment transport capacity of the subglacial water based upon the gradient of the hydraulic potential, channel size, and water flux (Table 1, Sect. 2.2; e.g., Walder and Fowler, 1994; Alley et al., 1997). The hydraulic model is based on the assumption that subglacial water flows along the hydraulic potential gradient, the weight of ice pressurizes water at the bed (Shreve, 1972), and the channel size varies over a substantially longer timescale compared to water discharge. This model includes characteristics of a Röthlisberger channel without explicitly describing properties such as creep closure and pressure melt of channel walls (Röthlisberger, 1972).

Table 1 Model variables.

The gradient of the hydraulic potential of a subglacial channel Ψ (at a certain location and time) can be determined with a known hydraulic diameter D_h (a function of channel size and shape) and water discharge Q_w. The gradient of the hydraulic potential can then be determined using the Darcy–Weisbach equation for fluid flow through a pipe:

$$\begin{array}{}\text{(1)}& \mathrm{\Psi }=s{f}_{\mathrm{r}}{\mathit{\rho }}_{\mathrm{w}}{\displaystyle \frac{Q_{\mathrm{w}}^{\mathrm{2}}}{D_{\mathrm{h}}^{\mathrm{5}}}},\end{array}$$

where the density of water is ρ_w, the Darcy–Weisbach friction factor is f_r, and the channel's cross-sectional geometry, which impacts water pressure, is accounted for by s (Hooke et al., 1990). We represent s as

$$\begin{array}{}\text{(2)}& s={\displaystyle \frac{\mathrm{2}(\mathit{\beta }-\sin \mathit{\beta }{)}^{\mathrm{2}}}{{\left(\frac{\mathit{\beta }}{\mathrm{2}}+\sin \frac{\mathit{\beta }}{\mathrm{2}}\right)}^{\mathrm{4}}}},\end{array}$$

where β is the central angle of the circular segment representing the channel edge. Smaller values of β result in broad channels and β=π results in a semicircular channel.

We assume that the hydraulic diameter D_h of the channel results from a characteristic water discharge $Q_{\mathrm{w}}^{*}$, which is evaluated by the source percentile of water discharge over a certain time period s_p and a response time of the channel size s_a, that remains consistent throughout the model run (Table 1; Delaney et al., 2019).

We sum the prescribed melt rate ${\dot{m}}_{\mathrm{w}}$ up the glacier to define Q_w, not considering englacial water storage. Percentile s_p over a response time period prior to the time step s_a is applied to Q_w to evaluate a characteristic water discharge $Q_{\mathrm{w}}^{*}$ that represents the size of the conduit (hours to days; see Gimbert et al., 2016; de Fleurian et al., 2018; Delaney et al., 2019; Nanni et al., 2020). The timescales, s_a, and characteristic water discharges (s_p and $Q_{\mathrm{w}}^{*}$) responsible for changes in subglacial conduit size are poorly constrained, yet their impact can be intuited. For instance, short-lived increases in water discharge due to an hour of precipitation will not greatly impact the hydraulic diameter of the subglacial channel, whereas prolonged melt would increase the hydraulic diameter.

With data of representative water discharge below the glacier $Q_{\mathrm{w}}^{*}$ and the static hydraulic pressure gradient Ψ^*, a representative hydraulic diameter D_h can be estimated. For a short period, such a D_h is assumed to be time-independent and is defined in Eq. (1):

$$\begin{array}{}\text{(3)}& D_{\mathrm{h}}={\left(s{f}_{\mathrm{r}}{\mathit{\rho }}_{\mathrm{w}}{{\displaystyle \frac{Q_{\mathrm{w}}^{*\mathrm{2}}}{\mathrm{\Psi }}}}^{*}\right)}^{\frac{\mathrm{1}}{\mathrm{5}}},\end{array}$$

where Ψ^* is a representative gradient of the hydraulic potential at overburden pressure, evaluated using the Shreve potential gradient:

$$\begin{array}{}\text{(4)}& {\mathrm{\Psi }}^{*}=\mathrm{\nabla }\left({\mathit{\rho }}_{\mathrm{i}}g\left(z_{\mathrm{s}}-z_{\mathrm{b}}\right)+{\mathit{\rho }}_{\mathrm{w}}g{z}_{\mathrm{b}}\right),\end{array}$$

where z_s and z_b are surface and bed elevations, respectively, ρ_i is the density of ice, and g is the gravitational acceleration constant.

With knowledge of D_h, we insert the instantaneous value of Q_w into Eq. (1) to evaluate the instantaneous gradient of the hydraulic potential Ψ. To prevent unreasonable water pressures when $Q_{\mathrm{w}}^{*}$ rapidly increases and D_h is small, the model limits the minimal hydraulic diameter to 0.3 m (Delaney et al., 2019).

2.2 Till layer model: bedrock erosion and sediment transport

SUGSET_2D simulates the evolution of a subglacial till layer, which we define as transportable sediment below the glacier produced through glacier erosion and fluvial sediment transport. Fluvial sediment transport, in supply- and transport-limited regimes, mobilizes and deposits sediment, removing or adding material from the till layer (Brinkerhoff et al., 2017; Delaney et al., 2019). Conversely, erosive processes such as abrasion and quarrying add material to the till layer. Note that we do not consider processes such as fluvial abrasion that appear to produce minimal sediment (Beaud et al., 2018b). To represent these processes, we implement the Exner equation (Fig. 2; Exner, 1920a, b; Paola and Voller, 2005), a mass conservation relationship, to solve for the till layer height given the erosive and fluvial conditions.

$$\begin{array}{}\text{(5)}& \underset{\mathrm{till}\mathrm{evolution}}{\underbrace{{\displaystyle \frac{\partial H}{\partial t}}}}=-\underset{\mathrm{sediment}\mathrm{transport}}{\underbrace{\mathrm{\nabla }\cdot Q_{\mathrm{s}}}}+\underset{\mathrm{bedrock}\mathrm{erosion}}{\underbrace{{\dot{m}}_{\mathrm{t}}}}\end{array}$$

H is till thickness and t is time (Table 1). The first term on the right side of the equation represents fluvial sediment transport processes, where ∇⋅Q_s represents sediment mobilization or deposition. The second term on the right side of the equation captures bedrock erosion processes, where ${\dot{m}}_{\mathrm{t}}$ is a bedrock erosion rate.

Figure 2 Illustration of a model cell (a), detailing the layers of bedrock, till, water, and ice. Characteristics of the subglacial channel are noted as a polygon but shown in one dimension for clarity in (b) with the Hooke angle parameterization with two different channel shapes for different values of β (Eqs. 2, 9, and 10).

We calculate sediment mobilization in both supply- and transport-limited conditions. Divergence of the sediment flux is evaluated by approximating ∇⋅Q_s with $\frac{{\stackrel{\mathrm{̃}}{Q}}_{\mathrm{s}}}{w}$ using a similar mobilization scheme as in Delaney et al. (2019)

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& {\stackrel{\mathrm{̃}}{Q}}_{\mathrm{s}}=\\ & \left\{\begin{array}{llll}\frac{Q_{\mathrm{sc}}-Q_{\mathrm{s}}}{l}& \mathrm{if}\frac{Q_{\mathrm{sc}}-Q_{\mathrm{s}}}{l}\le {\dot{m}}_{\mathrm{t}}w& \left(\mathrm{transport}\text{-}\mathrm{limited}\right)& \left(\mathrm{6}a\right)\\ \mathrm{0}& \mathrm{if}H=H_{\lim }\mathit{\&}\frac{Q_{\mathrm{sc}}-Q_{\mathrm{s}}}{l}\le \mathrm{0}& & \left(\mathrm{6}b\right)\\ \frac{Q_{\mathrm{sc}}-Q_{\mathrm{s}}}{l}\mathit{\sigma }\left(H\right)+{\dot{m}}_{\mathrm{t}}w(\mathrm{1}-\mathit{\sigma }(H\left)\right)& \mathrm{otherwise}& \left(\mathrm{supply}\text{-}\mathrm{limited}\right)& \left(\mathrm{6}c\right)\end{array}\right.\end{array}\end{array}$$

${\stackrel{\mathrm{̃}}{Q}}_{\mathrm{s}}$ is sediment mobilization across a width of the glacier bed w perpendicular to the water's flow direction. Note that w is not necessarily the channel width, but rather a representative width across the glacier bed over which sediment can be accessed by water flowing through the subglacial channel (Fig. 2). The channel width w_c is used to calculate the width over which to apply the sediment discharge capacity and is discussed below in Eq. (9), which converts the hydraulic diameter D_h to channel width. Q_sc is the sediment transport capacity or the maximum amount of sediment that could be transported under the given hydraulic conditions. l is a characteristic length scale for sediment mobilization, over which sediment mobilization adjusts to sediment transport conditions. σ is a sigmoidal function of H,

$$\begin{array}{}\text{(7)}& \mathit{\sigma }\left(H\right)={\left(\mathrm{1}+\exp \left(\mathrm{10}-\mathrm{5}{\displaystyle \frac{H}{\mathrm{\Delta }\mathit{\sigma }}}\right)\right)}^{-\mathrm{1}},\end{array}$$

which enables a smooth transition from transport- to supply-limited transport in Eq. (6c). If H, the till thickness, is greater than 3Δσ, then the impact on sediment mobilization is negligible and the system is in a transport-limited regime. As H approaches Δσ, then σ(H) is close to 0 and sediment transport is in a supply-limited regime since no significant sediment mobilization takes place.

The condition in Eq. (6a) represents the case in which bedrock erosion exceeds sediment mobilization, and thus sediment transport exists in a transport-limited regime. The condition in Eq. (6b) impedes mobilization or deposition, therefore transporting sediment to the next cell when a till thickness is equal to H_lim, the value of which is chosen to be on the order of the maximal change in till height over the model run (∼10 cm). The condition in Eq. (6b) prevents unbounded sediment accumulation, as the model does not include physical processes to limit sediment deposition, such as reduced channel size and increased water velocity in response to the infill of sediment (Perolo et al., 2018). The condition in Eq. (6c) allows sediment mobilization to transition between transport- and supply-limited regimes, limiting sediment mobilization to the sediment production term, ${\dot{m}}_{\mathrm{t}}$ (see below), when H is small and thus minimal sediment is available for transport. With the conditions in Eq. (6a–c), we can calculate sediment transport in transport- and supply-limited regimes and pass sediment through the system.

We calculate sediment transport capacity Q_sc using the total sediment transport relationship by Engelund and Hansen (1967),

$$\begin{array}{}\text{(8)}& Q_{\mathrm{sc}}={\displaystyle \frac{\mathrm{0.4}}{f_{\mathrm{r}}}}{\displaystyle \frac{\mathrm{1}}{D_{\mathrm{m}}{\left(\frac{{\mathit{\rho }}_{\mathrm{s}}}{{\mathit{\rho }}_{\mathrm{w}}}-\mathrm{1}\right)}^{\mathrm{2}}{g}^{\mathrm{2}}}}{\left({\displaystyle \frac{\mathit{\tau }}{{\mathit{\rho }}_{\mathrm{w}}}}\right)}^{\frac{\mathrm{5}}{\mathrm{2}}}{w}_{\mathrm{c}},\end{array}$$

where ρ_s (ρ_w) is the bulk density of the sediment (water), D_m is the mean sediment grain size, and τ represents the shear stress between the water and the channel bed.

The width of the channel floor w_c, required to evaluate the surface over which sediment transport may occur, is given by

$$\begin{array}{}\text{(9)}& w_{\mathrm{c}}=\mathrm{2}\sin {\displaystyle \frac{\mathit{\beta }}{\mathrm{2}}}\sqrt{{\displaystyle \frac{\mathrm{2}S}{\mathit{\beta }-\sin \mathit{\beta }}}},\end{array}$$

where, again, β is the Hooke angle controlling channel morphology (Sect. 2.1), and S is the cross-sectional area of the channel given by

$$\begin{array}{}\text{(10)}& S={\displaystyle \frac{D_{\mathrm{h}}^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{{\left(\frac{\mathit{\beta }}{\mathrm{2}}+\sin \frac{\mathit{\beta }}{\mathrm{2}}\right)}^{\mathrm{2}}}{\mathit{\beta }-\sin \mathit{\beta }}}.\end{array}$$

Here, hydraulic diameter D_h is evaluated from Eq. (3).

We also determine the shear stress between water flowing through the channel and the sediment below in Eq. (8) through the Darcy–Weisbach relationship,

$$\begin{array}{}\text{(11)}& \mathit{\tau }={\displaystyle \frac{\mathrm{1}}{\mathrm{8}}}{f}_{\mathrm{r}}{\mathit{\rho }}_{\mathrm{w}}{v}^{\mathrm{2}},\end{array}$$

where $v=\frac{Q_{\mathrm{w}}}{S}$ is the water velocity. Water discharge Q_w is calculated by the water flowing above a position in the glacier, and S, the cross-sectional area, is evaluated in Eq. (10). Other sediment transport relationships using shear stress could be exchanged by the model operator (e.g., Meyer-Peter and Müller, 1948). We have chosen the Engelund and Hansen (1967) formulation due to the representation of both suspended and bedload transport.

We assume that till armors the bed against erosion (e.g., Alley et al., 2003; Brinkerhoff et al., 2017; Delaney et al., 2019). In response, the source term, ${\dot{m}}_{\mathrm{t}}$, is described as

$$\begin{array}{}\text{(12)}& {\dot{m}}_{\mathrm{t}}=\dot{e}\left(\mathrm{1}-{\displaystyle \frac{H}{H_{\max }}}\right),\end{array}$$

where H_max is a till height beyond which no further erosion, $\dot{e}$, may occur.

We use an empirical relationship with sliding velocity u_b to describe bedrock erosion,

$$\begin{array}{}\text{(13)}& \dot{e}=k_{\mathrm{g}}{u}_{\mathrm{b}}^{l_{\mathrm{er}}},\end{array}$$

where k_g is an erodibility constant and l_er is an exponent, which varies between 0.66 and 3 (Herman et al., 2021). The sliding velocity, u_b, is assumed to be related to basal shear stress (τ_b;  Weertman, 1957) given the following relationship,

$$\begin{array}{}\text{(14)}& u_{\mathrm{b}}=B{\mathit{\tau }}_{\mathrm{b}}^{\mathrm{m}},\end{array}$$

where B is a constant and we assume the exponent m is equal to 1.

We assume that τ_b is equal to driving stress (Cuffey and Paterson, 2010),

$$\begin{array}{}\text{(15)}& {\mathit{\tau }}_{\mathrm{b}}={\mathit{\rho }}_{\mathrm{i}}gh\sin \left(\mathit{\alpha }\right),\end{array}$$

where ρ_i is the density of ice, h is the glacier thickness, and α is the surface slope of the glacier.

Note that alternative parameterizations of erosion or basal sliding can easily be exchanged for ${\dot{m}}_{\mathrm{t}}$.

2.3 Spatial and temporal discretization as well as parameters

We describe the water and sediment routing and numerical implementation of the equations presented above.

2.3.1 Water and sediment routing and implementation

A routing scheme is implemented to (1) evaluate the hydraulic potential and thus the direction of the water flow and (2) transport sediment and water across the glacier bed to where it is expelled or deposited.

Figure 3 Example of model parameters and variables for the snapshot of the Griesgletscher case in Sect. 3.2. (a) Glacier flotation fraction. (b) Channel cross-sectional area S with (d) distributed water discharge, (c) the number of receiver cells r_t for a given cell, and (e) the water velocity.

Figure 4 Routing scheme on the grid. Solid lines represent cell boundaries, blue squares are cell centers, and red squares are cell edges. ϕ, the hydraulic potential, decreases in the direction of arrows so that water and sediment generally flow left to right and top to bottom. Edge length (λ) and cell area (δ) are shown. Cell numbers refer to identification in the stack (s_t). Select cells denote the weight of donors $w_{\mathrm{d},i,j}$, number of donors n_d, donor cells d_n, number of receivers n_r, and receiver cells r_s. Variables and their respective locations on the grid are shown. Some red and blue squares have been removed in some cells for clarity.

To evaluate the hydraulic potential and thus the direction of the water and sediment flow across the glacier bed, we use a two-dimensional routing scheme (Quinn et al., 1991) to establish flow routing based on the steepest hydraulic potential. We implement this scheme in a similar fashion as Bovy et al. (2016), but on a regular grid with square cells, extending in x and y directions. Water and sediment fluxes can pass to the four surrounding cells sharing an edge as a result of the x and y components of the hydraulic gradient at a given point in time. This routing algorithm returns a stack (s_t; Table 3), which is a vector that contains information about the order of cells to perform the calculations, along with the number of cells flowing into a cell (donors; n_d), the number of cells to which a single cell contributes (receivers; n_r), and the weight or the percentage of hydraulic potential and water (or sediment) discharge directed from one cell to another (w_d or w_r), as determined by the hydraulic potential gradient between cells (Fig. 4).

We define the hydraulic potential, upon which the routing scheme is evaluated, at a cell i in the grid, ϕ_i, based upon the elevation of the glacier bed plus the ice thickness, following Shreve (1972):

$$\begin{array}{}\text{(16)}& {\mathit{\varphi }}_i=f_{\mathrm{f}}{\mathit{\rho }}_{\mathrm{i}}g\left(z_{\mathrm{s},i}-z_{\mathrm{b},i}\right)+{\mathit{\rho }}_{\mathrm{w}}g{z}_{\mathrm{b},i},\end{array}$$

where f_f is the flotation fraction across the glacier, z_s is the glacier surface, and z_b is the glacier bed.

For the initial time step, the hydraulic potential ϕ is evaluated under the condition that f_f=1. After the first time step, we assume that the flotation fraction will vary in response to changing hydraulic conditions, such as diurnal or seasonal water input (e.g., Iken and Bindschadler, 1986). In turn, to establish an average flotation fraction, f_f, across the glacier bed for Eq. (16), we use

$$\begin{array}{}\text{(17)}& f_{\mathrm{f}}=\mathrm{mean}\left({\displaystyle \frac{{\mathit{\varphi }}_{\mathrm{0},i}}{{\mathit{\rho }}_{\mathrm{i}}g\left(z_{\mathrm{s},i}-z_{\mathrm{b},i}\right)+{\mathit{\rho }}_{\mathrm{w}}g{z}_{\mathrm{b},i}}}\right),\end{array}$$

where the denominator represents the hydraulic potential at the overburden pressure of the glacier (f_f=1 in Eq. 16), and i represents a cell in the grid.

ϕ₀ represents the hydraulic potential evaluated from summing the gradient of the hydraulic potential, Ψ, in Eq. (1) up-glacier. ϕ₀ at each cell i is evaluated as

$$\begin{array}{}\text{(18)}& {\mathit{\varphi }}_{\mathrm{0},i}={\mathrm{\Psi }}_{i,j}\cdot \mathit{\lambda }+\sum _{j=\mathrm{1}}^{n_{\mathrm{r}}}\left({\mathit{\varphi }}_{\mathrm{0},j}\cdot w_{\mathrm{r},i,j}\right).\end{array}$$

Here, by evaluating Eq. (1), Ψ_i,j is established for receiver cell j of i, λ is the edge length of a cell on a regular grid, n_r is the number of receivers that the cell i has, and $w_{\mathrm{r},i,j}$ is the proportion of the hydraulic potential fed by the upstream cell j to cell i. The operation is executed on a cell-by-cell basis, beginning with cells that have no receivers, such as those near the glacier terminus, and moving up the glacier using the inverted stack in s_t (Fig. 4).

Using the routing scheme above, we evaluate the water discharge from cell i, Q_w,i, from melt upstream as

$$\begin{array}{}\text{(19)}& Q_{\mathrm{w},i}=\sum _{j=\mathrm{1}}^{n_{\mathrm{d}}}{Q}_{\mathrm{w},j}\cdot w_{\mathrm{d},i,j}+{\dot{m}}_{\mathrm{w},i}\cdot \mathit{\delta },\end{array}$$

where n_d is the number of donor cells for cell i, $w_{\mathrm{d},i,j}$ is the percentage of water flow from cell j to cell i, and ${\dot{m}}_{\mathrm{w},i}$ is a prescribed meltwater source term in cell i. The operation begins with cells that have no donors (for instance at the top of the glacier) so that water accumulates down the glacier (Fig. 4). The amount of sediment leaving a cell i, Q_s,i, is the flux into the cell plus the sediment mobilized in the cell, which is defined as

$$\begin{array}{}\text{(20)}& Q_{\mathrm{s},i}=\sum _{j=\mathrm{1}}^{n_{\mathrm{d}}}{Q}_{\mathrm{s},j}\cdot w_{\mathrm{d},i,j}+{\stackrel{\mathrm{̃}}{Q}}_{\mathrm{s},i}\cdot \mathit{\lambda }.\end{array}$$

The first term is the flux of sediment into the cell i from donor cells j. The second term is sediment mobilization, ${\stackrel{\mathrm{̃}}{Q}}_{\mathrm{s},i}$, in cell i, which is computed by implementing Eq. (6a–c) as follows.

$$\begin{array}{}\text{(21)}& \begin{array}{rl}& {\stackrel{\mathrm{̃}}{Q}}_{\mathrm{s},i}=\\ & \left\{\begin{array}{llc}\sum _{j=\mathrm{1}}^{n_{\mathrm{d}}}\left(\frac{Q_{\mathrm{sc},j}-Q_{\mathrm{s},j}}{l}\cdot w_{\mathrm{d},i,j}\right)& \mathrm{if}\sum _{j=\mathrm{1}}^{n_{\mathrm{d}}}\left(\frac{Q_{\mathrm{sc},i,j}-Q_{\mathrm{s},i,j}}{l}\right)\cdot w_{\mathrm{d},i,j}\le \frac{{\dot{m}}_{\mathrm{t},i}\mathit{\lambda }}{n_{\mathrm{d}}}& \left(\mathrm{21}\mathrm{a}\right)\\ \mathrm{0}& \mathrm{if}{H}_j=H_{\lim }\mathit{\&}\frac{Q_{\mathrm{sc},i,j}-Q_{\mathrm{s},i,j}}{l}\le \mathrm{0}& \left(\mathrm{21}\mathrm{b}\right)\\ \frac{{\dot{m}}_{\mathrm{t},i}\mathit{\lambda }}{n_{\mathrm{d}}}(\mathrm{1}-\mathit{\sigma }(H\left)\right)+\sum _{j=\mathrm{1}}^{n_{\mathrm{d}}}\left(\frac{Q_{\mathrm{sc},j}-Q_{\mathrm{s},j}}{l}\right)\cdot \mathit{\sigma }\left(H\right)\cdot w_{\mathrm{d},i,j}& \mathrm{otherwise}& \left(\mathrm{21}\mathrm{c}\right)\end{array}\right.\end{array}\end{array}$$

Q_sc,j is the sediment transport capacity from cell j flowing to i, Q_s,j is sediment discharge entering from cell j to cell i, l is a response length scale, and λ is edge length.

We evaluate the change in till height at a cell by implementing Eq. (5) as

$$\begin{array}{}\text{(22)}& {\displaystyle \frac{\mathrm{d}{H}_i}{\mathrm{d}t}}={\displaystyle \frac{-Q_{\mathrm{s},i}+\sum _{j=\mathrm{1}}^{n_{\mathrm{d}}}{Q}_{\mathrm{s},j}\cdot w_{\mathrm{d},i,j}}{\mathit{\delta }}}+{\dot{m}}_{\mathrm{t},i},\end{array}$$

where δ is the cell area (Fig. 4). The term Q_s,i is the amount of sediment leaving the cell from Eq. (20), and the term $\sum _{j=\mathrm{1}}^{n_{\mathrm{d}}}{Q}_{\mathrm{s},j}\cdot w_{\mathrm{d},i,j}$ is the sediment flux entering the cell from the donors.

2.3.2 Numerics and parameters

Spatial discretization on the regular grid must be substantially smaller than characteristic length scale, l, in Eqs. (6a–c) and (21). We then solve Eq. (22) to establish till height, H, for given initial and boundary conditions in response to till production, ${\dot{m}}_{\mathrm{t}}$, and the divergence of the sediment discharge, Q_s, using an explicit time integration scheme.

To discretize the problem in time, the model implements the VCABM solver (Hairer et al., 1992; Radhakrishnan and Hindmarsh, 1993) from the package DifferentialEquations.jl (Rackauckas and Nie, 2017) to evolve till layer height, H. This solver implements an adaptive time step and uses a linear Adams–Moulton multistep method that is well-suited to non-stiff problems, which is optimal because of the rapid fluctuations in sediment transport that can occur. We impose a maximum time step of 6 h to ensure that the model captures the response to diurnal variations in melt input. In practice, the solver commonly uses a time step of roughly 20 min, which varies depending on sediment transport conditions and solver tolerance. Longer time steps occur over periods with minimal glacier melt, and thus sediment transport ceases (i.e., winter months). Table 3 presents the numerical parameters used.

We execute the routing scheme based upon hydraulic conditions to the nearest 6 min to improve stability and fill closed basins in the hydraulic potential to maintain continuous sediment transport through the domain. Smaller solving tolerances increase the computational time due to (1) the increased accuracy of the solution and (2) the reassessment of flow fractions between the adjacent cells, which results in different routing configurations as the model converges.

We impose boundary conditions on the edge cells so no sediment or water enters the domain. At outlet cells, water discharge leaves the domain, as does a flux of sediment, based on sediment transport conditions. In other applications, boundary conditions could also be set to represent processes such as hillslope erosion or glacial lakes, which route sediment or water to the subglacial environment (e.g., Andersen et al., 2015). At the outlet cells, we assume that the hydraulic potential has no ice overburden pressure.

Evolving Eq. (5) requires an initial till height, H₀, chosen by the model user. This initial till height represents material from bedrock erosion created prior to the model initialization. We apply a “spin-up” procedure to create a reasonable relationship between the amount of fluvial sediment transport and bedrock erosion.

New versions of the code are tested against reference cases to ensure consistency. Additionally, in each test, we ensure mass conservation by verifying that the amount of sediment leaving the system through fluvial transport is consistent with the till height change and erosion occurring under the simulated glacier.

3 Model application

We use two case studies to highlight model viability under increasingly complex situations. First, we apply the model to a synthetic alpine glacier topography with a synthetic hydrologic forcing, based on the Subglacial Hydrology Model Intercomparison Project (SHMIP; de Fleurian et al., 2018), to illustrate the model's performance in a simplistic scenario. We then apply the model to the topography, sediment, and water discharge at Griesgletscher in the Swiss Alps. We demonstrate its proficiency by comparing the calculated sediment transport output with measured data (Delaney et al., 2018a). We also identify important factors controlling subglacial sediment discharge in these simulations.

3.1 Synthetic alpine case

3.1.1 Experiment design

We run simulations using a synthetic alpine glacier geometry, along with seasonally and diurnally varying hydrological forcing from the SHMIP experiments (de Fleurian et al., 2018). The domain is 6000 m on one axis and 1080 m on the other (Fig. 7). The resulting geometry approximates the Bench Glacier in Alaska. The U-shaped bed and variable ice thickness mean that variable hydrologic gradients will occur perpendicular to the flow, and thus water and sediment are routed across multiple cells.

To represent a hydrological environment that varies both seasonally and diurnally, we implement a simple spatially distributed melt model, as in SHMIP (de Fleurian et al., 2018):

$$\begin{array}{}\text{(23)}& {\dot{m}}_{\mathrm{w}}\left(z_{\mathrm{s}}\right)=\left\{\begin{array}{ll}{M}_{\mathrm{f}}T\left(z_{\mathrm{s}}\right)+{\dot{m}}_{\mathrm{b}}& \mathrm{if}T\left(z_{\mathrm{s}}\right)>\mathrm{0}\\ {\dot{m}}_{\mathrm{b}}& \mathrm{if}T\left(z_{\mathrm{s}}\right)\le \mathrm{0}\end{array}\right.,\end{array}$$

where M_f=0.01 m ^∘C⁻¹ d⁻¹ is a melt factor, and ${\dot{m}}_{\mathrm{b}}$ is the basal melt rate. T(z_s) is air temperature T (^∘C) at elevation z_s, defined as

$$\begin{array}{}\text{(24)}& \begin{array}{rl}T\left(z_{\mathrm{s}}\right)& =\left(-A_{\mathrm{a}}\cos \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }t}{s_{\mathrm{year}}}}\right)+A_{\mathrm{d}}\cos \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }t}{s_{\mathrm{day}}}}\right)+\mathrm{\Delta }T-\mathrm{5}\right)\\ & \cdot \left(\mathrm{1}+z_{\mathrm{s}}{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}z}}\right),\end{array}\end{array}$$

where A_a and A_d are the annual and diurnal amplitudes in temperature, respectively, ΔT is a temperature offset that is adjusted to control the meltwater input, s_day represents the number of seconds in 1 d, s_year is the number of seconds in a year, and $\frac{\mathrm{d}T}{\mathrm{d}z}=-\mathrm{0.0075}$ ^∘C m⁻¹ is the air temperature lapse rate (de Fleurian et al., 2018). In this case, we route water directly to the subglacial system at the location where the melt occurs, omitting moulins or crevasses that concentrate meltwater delivery to the bed, for instance.

We run the model for 10 years with a steady climate; we then apply a linear temperature increase of 0.5 ^∘C a⁻¹ for 10 years followed by 10 years of steady temperature at the maximal ΔT. We implement this dramatic warming to capture the model's ability to represent different climatic conditions. The model is initiated with 5 cm of till across the bed. To spin up the model, we apply the initial year of hydrological forcing for 5 years to limit computational time. In other applications, the spin-up could be maintained until the annual change in till height was well below commonly accepted glacier erosion rates (Hallet et al., 1996).

3.1.2 Model outputs and findings

The simulations show that over seasonal timescales, sediment discharge increases at the onset of melt and decreases shortly thereafter, but prior to the maximum amount of water discharge that occurs at the peak of each melt season (Fig. 5). Maximum and average quantities of daily sediment discharge decrease until the very end of the melt season, when sediment discharge increases very slightly again (Figs. 5 and 6b, d, f). This slight increase occurs when water stops flowing during the night, allowing sediment from bedrock erosion to briefly accumulate in the channels from bedrock erosion. Increased sediment discharge produced by the model at the beginning of the melt season results from greater sediment availability following the growth of the till layer over the winter months when the small amount of melt prevents substantial transport of sediment. In fact, similar increases in sediment discharge have been observed in alpine glaciers at the onset of melt each season (Figs. 5b and 6b, d, f) (Willis et al., 1996; Swift et al., 2005; Riihimaki et al., 2005; Delaney et al., 2018b) and are reproduced in the one-dimensional version of the model in Delaney et al. (2019).

Figure 5 Model output from a synthetic alpine topography and forcing over a 30-year run with diurnal and seasonal variations in melt input as daily mean values. The gray box represents a time period of increasing glacier melt. (a) Daily averaged seasonally varying water discharge (Q_w) increases from year 10 to 20, while till height (H) decreases throughout the model run, with seasonal increases in the absence of glacier melt. (b) Daily averaged sediment discharge in brown shows strong seasonal variability. Annual sediment discharge (green) increases with increasing melt, with the highest sediment discharge occurring in year 19 when glacier melt is greatest. During stable climate temperatures before and after the increase in temperature, annual sediment discharge generally decreases. However, following the melt, Q_s stabilizes at a higher level due to the increased area over which sediment is transported.

Figure 6 Annual response to different till production patterns across the synthetic glacier case studies. (a, b) Conventional model setup in which sediment is produced year-round, ORIGINAL. (c, d) Setup equivalent to the previous setup, except sediment is produced only in summer months when water is present at the glacier bed, SEASON. Note that till height remains constant on the edges of the plot over the winter months. (e, f) Steady erosion of 2 mm a⁻¹ across the entire synthetic glacier, with no spatial or temporal variability in sediment production, CONST. Data are plotted at a 6 h interval so that daily maximums and minimums are visible.

Over the course of the simulation, the mean till height across the glacier generally decreases throughout the model run, although there are small increases in till height during the winter months without sediment transport and larger decreases in till height during the time periods with substantial melt (Figs. 5a and 6a, c, e). Exhaustion of sediment is evident in the middle of the glacier where much of the water flows as a result of the spin-up procedure prior to the model initialization (Fig. 7e and f). Note that the decreasing till height through the model run results from sediment mobilization on the margins of the glacier, where increased water flow occurs more often in a warmer climate. During the climate warming from years 10–20, sediment discharge from the glacier increases due to greater melt and water discharge on the upper reaches of the glacier. This results in increased sediment transport at higher elevations on the glacier, where sediment could persist in a cooler climate (Fig. 7e and f). Following the stabilization of the climate at year 20, sediment discharge remains elevated compared to the cooler climate because sediment transport occurs over a larger region of the glacier bed (Fig. 5b).

Figure 7 Spatial view of synthetic case (Sect. 3.1). Spatially distributed (a) ice thickness (h), (b) sediment discharge (Q_s), (c) water discharge (Q_w), (d) glacier erosion rate ($\dot{e}$), (e) till layer height (H) at year 10 prior to warming, and (f) till layer height (H) at year 29 near the end of the model run. Spatial differences in the distribution of water and sediment discharge in panels (e) and (f) result from the depletion of subglacial till beneath the glacier. We have included an animation of this figure in the Video supplement.

In the model, bedrock erosion relies only on driving stress and till thickness. Sliding and bedrock erosion did not vary seasonally with increased subglacial water discharge (Fig. 5a). This causes sediment to accumulate during the winter months, which subsequently provides ample material for transport when melting increases in the spring. To test the effects of spatially variable erosion and the role of hydrology, we present two additional cases to supplement the synthetic alpine glacier case above, named ORIGINAL. An additional synthetic case, SEASON, simulates bedrock erosion by only allowing sliding, and thus erosion, during the summer months (e.g., Iken and Bindschadler, 1986; Herman et al., 2011); the same erosional relationship is applied as that in Sect. 3.1 (Eqs. 13 and 12). In the SEASON case, however, erosion occurs only when the amount of water input substantially exceeds the background basal melt input rate that is present in the winter. This case captures the seasonal variations in bedrock erosion (Ugelvig et al., 2018). In another additional case, CONST, bedrock erosion remains constant over the entirety of the glacier at a rate of 2 mm a⁻¹, independent of the spatially varying glacier sliding velocity of the other cases (Fig. 7).

The ORIGINAL case discharges over 11 620 m³ of sediment per year, while the SEASON case discharges only 60 % of that value due to the absence of bedrock erosion during the winter months. The CONST case discharged 7320 m³ of sediment over the year or 63 % of the ORIGINAL case. The quantity of sediment discharge in the CONST case results in a catchment-scaled height change in the till layer of approximately 1.1 mm a⁻¹ due to decreased erosion efficiency with till height instead of the prescribed bedrock erosion rate of 2 mm a⁻¹ (Eq. 12). Additionally, the spatial disparity of where sediment is produced at the glacier bed compared to the location of sediment transport further reduces the catchment-scaled height change (Fig. 7e and f).

In each of the three synthetic cases, sediment discharge increases at the onset of melt and substantially decreases by the end of the melt season due to sediment exhaustion (Fig. 7). In ORIGINAL (Fig. 6a and b), greater sediment discharge occurs compared to the alternate cases (SEASON and CONST). The increased sediment discharge in ORIGINAL results from (1) the winter periods without melt over which bedrock erosion adds more sediment to the layer in the absence of sediment transport and (2) bedrock erosion that occurs low on the glacier where much of the sediment transport occurs (Fig. 7d). By contrast, in the CONST case, a steady amount of erosion occurs across the entire glacier bed. The peak sediment discharge in CONST (Fig. 6e and f) occurs slightly earlier in the season compared to the ORIGINAL and SEASON cases due to the increased amounts of sediment below the glacier's lower portions.

3.2 Griesgletscher

3.2.1 Experiment design

We also run simulations of Griesgletscher in the Swiss Alps using topographic data from Delaney et al. (2019). Here, we run the model from 2009–2017, with the modeled water discharge time series from Delaney et al. (2018a). Subglacial sediment discharge from the glacier is determined over four different time periods (2011–2013, 2013–2014, 2014–2015, 2015–2016) by differencing the bathymetry maps collected through this period and considering proglacial erosion quantities (Delaney et al., 2018a, 2019). To estimate surface melt across the glacier with respect to elevation, we use

$$\begin{array}{}\text{(25)}& {\dot{m}}_{\mathrm{w}}(x,y)={\dot{b}}^{\mathrm{0}}+\mathit{\gamma }\left(z_{\mathrm{s}}(x,y)-z_{\mathrm{s}}^{\mathrm{0}}\right).\end{array}$$

Here, γ is the mass balance gradient, $z_{\mathrm{s}}^{\mathrm{0}}$ represents the glacier's lowest elevation, and ${\dot{b}}^{\mathrm{0}}$ represents the melt rate at the glacier's lowest extent. ${\dot{b}}^{\mathrm{0}}$ was evaluated numerically at each water discharge value using the hypsometry of the glacier.

We apply a parameter search over a range of values of sediment grain size (D_m, representing a primary control on fluvial transport of subglacial sediment), sliding rate factor (B, representing a control on bedrock erosion), and the initial till height condition (H₀, representing the effects of existing quantities of sediment below the glacier). A total of 100 simulations were run with randomly selected parameters from a uniform distribution. No spin-up was applied in this case to establish an initial condition because of the wide range of H₀ values explored.

The wall time for a single model run averaged 8.9 h, and each run for a parameter set was executed on a single CPU. Instead of applying the mean flotation fraction across the glacier, as done in the synthetic cases, the maximum value was applied with an upper limit of 1.

We consider model outputs resulting in a perfect rank correlation across the four data collection periods and have errors less than the 131 000 m³ of sediment that was expelled from the glacier over this period (Delaney et al., 2018a, 2019). For the example presented below, we show the simulation with the lowest absolute error between the model output and the sediment transport data.

3.2.2 Model outputs and findings

The parameter search yields an optimum grain size parameter D_m of 2 cm, sliding parameter B of $\mathrm{2.05}\times {\mathrm{10}}^{-\mathrm{11}}$ MPa m s⁻¹, and initial till height H₀ of 2.5 mm (see red stars, Fig. 8a–c). The model's ability to reproduce the quantities of sediment in the validation data largely depends on the grain size parameter, D_m, shown by Fig. 8a. Compared to D_m, the sliding parameters (B) and initial condition parameters (H₀) have a reduced influence in representing the measured data, given that similar values of B and H₀ can produce largely different model outputs in the context of D_m (Fig. 8a–c).

Figure 8 Results of the parameter search (a–c), the frequency of parameter values that produced a rank correlation of 1 (d–f), and average sediment flux from the model run amongst the parameter combinations over the time periods (g) in the Griesgletscher case. Red stars represent the optimum parameter combination with an absolute error of roughly 62 600 m³. Blue lines represent all model outputs, while the gray line represents the optimum parameter combination.

The optimized parameter combination, along with other parameter combinations, reproduces the interannual variability in sediment discharge from the Griesgletscher (Fig. 8g). The absolute error between the optimum model run and the measured data is roughly 62 600 m³ of sediment. The error is slightly less than half of the 131 300 m³ total sediment discharged from the Griesgletscher over this time period from 2011 to 2016 (Delaney et al., 2018a). In fact, the model runs capture the third period from late 2014 to late 2015 well. However, the runs systematically overestimate the second and fourth periods and generally underestimate the high-sediment-discharge period from late 2011 until late 2013 (Fig. 8g).

Figure 9 Time series of model output from optimum parameter combinations. (a) Daily averaged water discharge, an input modeled for Griesgletscher, Switzerland, in Delaney et al. (2018b), and sediment discharge, the output of the model, from Griesgletscher beginning in 2010. (b) Daily averaged sediment transport capacity and average till height. Note that sediment discharge capacity is roughly 1 order of magnitude larger than sediment transport discharge. Additionally, the increasing trend in till height, H, through this model run shows that sediment is produced at a greater rate than it is transported from the glacier bed.

The best-performing model run shows strong temporal variability in sediment discharge (Fig. 9a). Some of the peaks in sediment discharge occur during the short periodic increases in water discharge. Yet the greatest sediment discharge values do not necessarily occur at the highest water discharge values (Figs. 9a and 10a). Despite the strong dependence on grain size and fluvial transport of sediment in the parameter search (Fig. 8a), the modeled sediment transport capacity Q_sc still remains roughly an order of magnitude higher than the calculated sediment discharge Q_s (Fig. 9a and b). The steep section of the glacier (Fig. 11c) experiences sediment depletion over the model run, as do several patches of the glacier bed near the overdeepening and high on the glacier (lower left of panels in Fig. 11c and d). On some parts of the upper glacier, the calculated bedrock erosion grows the till layer beyond the initial condition in the absence of substantial sediment transport.

Table 2 Physical model parameters and constants.

Figure 10 Model outputs of sediment discharge from the glacier compared to water discharge (a) and sediment transport capacity (b).

The range of values for B, the sliding parameter, in the parameter search results in mean sliding velocities across the glacier between 14 and 70 m a⁻¹ (Eq. 14). The optimum run in the parameter search results in an average sliding velocity of 39 m a⁻¹. We note that smaller sliding velocities could result in equivalent amounts of erosion if the parameters k_g and l_er in Eq. (13) are increased. The model reveals that the relatively large velocities needed to produce adequate amounts of erosion occur in part because sediment production in the model is limited to the narrow patches of the glacier bed where minimal till persists and bedrock erosion may occur. As a result, the model requires more sliding to produce the equivalent amount of sediment with more till at the bed, even though the sliding and erosion parameters applied here are within a reasonable range (Fig. 11c and d).

4 Model limitations

The lack of observations of the spatial distribution of subglacial sediment makes selecting an initial value of H, the height of the dill layer, difficult. The slow rate of basal erosion suggests that an equilibrium between fluvial sediment transport and bedrock erosion will likely take centuries to attain if such an equilibrium could even exist in light of the variability in climatic, and thus glacier, conditions. Should an equilibrium eventually be present (e.g., Herman et al., 2018; Delaney and Adhikari, 2020), it is probably beyond a feasible computational time for this model, given its processing speeds.

Figure 11 Spatial view of characteristics from the Griesgletscher model run (see Fig. 1 images). (a) Subglacial sediment transport is concentrated in a narrow part of the bed, (b) and water discharge is also highly variable across the bed. Till layer height changes substantially from the beginning of the model run (c) to the end of the model run (d). We identify the overdeepening near the glacier terminus as well as a steep section connecting the upper and lower glacier. Over this time, till exhaustion in regions of high water flow is visible, while other regions of the glacier bed experiencing sediment deposition and till growth from glacier erosion can be identified. We have included an animation of this figure in the Video supplement.

SUGSET_2D also contains 20 parameters (Tables 2 and 3). In the available literature, some of these parameters have been partially constrained using inverse methods (Brinkerhoff et al., 2016) as well as detailed modeling and measurements (e.g., Chen et al., 2018; Covington et al., 2020; Pohle et al., 2022). However, many are poorly constrained.

Table 3 Numerical model parameters.

For instance, we limit the thickness at which the till layer must stop accumulating sediment (Eq. 6b, H_lim) due to the changes in the hydraulic potential caused by channel infill of sediment. We assume that this value is on the order of tens of centimeters (Table 2) based upon available observations of sediment deposition and glacier uplift (Perolo et al., 2018). While the impact of a till layer on bedrock abrasion remains uncertain, we expect that sediment of a certain thickness will armor the bed, preventing erosion (Alley et al., 2003). In turn, we limit erosion with till thickness to a threshold (5 cm) of the same order as H_lim to improve computational time. Due to the difficulty of making direct observations at glacier beds, only one study, to our knowledge, has quantified till thickness at a single point below a glacier (Truffer et al., 2000). The initial till height, H₀, in the model must therefore be chosen thoughtfully because the system will remain impacted by this boundary condition throughout the model run. Furthermore, the model does not consider the interactions between fluvial sediment transport, debris concentrations in subglacial ice, and bedrock erosion, which may be important for subglacial sediment transport (e.g., Ugelvig et al., 2018).

The routing method we use assumes that the water flow direction is in response to the Shreve potential (Sect. 2.3.1). Therefore, it does not explicitly simulate the evolution of efficient and inefficient subglacial drainage systems over the course of the season or the inheritance of existing subglacial canals or channels (Fig. 3; e.g., Werder et al., 2013; Zechmann et al., 2021). In addition, a response time of the subglacial channel is chosen prior to simulations; this could be compared to a more sophisticated, but computationally more expensive, representation of processes in an R-channel model (e.g., Röthlisberger, 1972). Lastly, we have chosen a single friction factor f_r for the entirety of the run. This factor can vary in time (Pohle et al., 2022) and can be impacted by other factors such as sediment grain size or bedrock along the channel bed.

5 Implications

Results of the one-dimensional model (SUGSET; Delaney et al., 2019) and the two-dimensional model, SUGSET_2D, highlight the importance of simulating the spatial heterogeneities in bedrock erosion, sediment availability, and sediment transport capacity. In the one-dimensional version of SUGSET, water can access sediment from the till layer across the entire glacier width, perpendicular to the glacier flow line. In SUGSET_2D, however, sediment access and transport are not averaged over the glacier width. Rather, by considering the spatial distribution in water discharge and sediment availability laterally below a glacier, the model evaluates where heterogeneities may persist and how these heterogeneities will impact subglacial sediment dynamics (Figs. 7 and 11).

In SUGSET_2D, large diurnal increases in sediment discharge occur near peak daily melt because the area of flowing water expands under the glacier (Fig. 6b, d, and f). As a result, increased sediment transport can occur in regions of the glacier bed with substantial sediment when hydraulic conditions permit; that patch of bed is abandoned when water is routed to another part of the glacier bed (Video supplement). This allows sediment to be stored in these regions of the bed until the hydraulic conditions return and renew the increase in sediment transport. Such processes cannot be represented in a one-dimensional model, wherein the entire width of the glacier evolves together (Figs. 5b and 6b, d, f; Video supplement).

When we compare the model runs across the space of three parameters, D_m, B, and H₀, to sediment discharge data from Griesgletscher in the Swiss Alps (Sect. 3.2) we find that the model has a limited ability to capture the large sediment discharge from the first time period (2011–2013). The reduced sediment discharge in the second and fourth time periods, 2013–2014 and 2015–2016, indicates that the model does not adequately represent the processes that are responsible for the measured increases in sediment transport in these time periods (Fig. 8). Such processes may include variable flow routing following channel shape (Eqs. 1 and 2) or flotation fraction (Eq. 16). Additional processes may also be omitted due to model inputs. For instance, the evolving surface topography, not considered here, may cause alternative flow paths below the glacier (Fischer et al., 2005) and exposes new patches of the glacier bed to sediment transport or the relocation of channels (e.g., Zechmann et al., 2021). Furthermore, glacier sliding remains constant over the model run, so the results do not explicitly account for seasonal or interannual variability in bedrock erosion (e.g., Herman et al., 2015; Ugelvig et al., 2018); however, temporal variations in bedrock erosion are calculated through changing till thickness (Eq. 12 and Fig. 9).

Model performance at Griesgletscher depends greatly on sediment grain size compared to other parameters such as the initial till condition or bedrock erosion (Fig. 8c and d). Grain size is a strong control on sediment discharge in SUGSET_2D because it modulates sediment mobilization in patches of the bed only occasionally accessed by subglacial flow during the melt season – after sediment has been largely evacuated from the main channel (Fig. 11). This process cannot be fully represented in a one-dimensional model. This is especially so because the main flow paths under the glacier can be evacuated of sediment (Fig. 11c and d). Thus, these flow paths contribute to the catchment's sediment discharge only through the new production of sediment through erosion (Eq. 12). However, the dependence on grain size suggests that the connectivity between subglacial channels and distal sediment patches is a strong control on sediment discharge from the subglacial system. While this seems to be an important process on the relatively small and shallow Griesgletscher alpine glacier (Delaney et al., 2018a), sediment discharge on other, potentially steeper, glaciers may respond more strongly to processes such as bedrock erosion (Herman et al., 2015). The connectivity between the main channels and distal sources of sediment occurs through the fluvial transport capacity of sediments here, but this connectivity may also be influenced through processes not considered in the model, such as till deformation (e.g., Damsgaard et al., 2020) or sediment sorting (e.g., Bacchi et al., 2014).

Lastly, the model demonstrates the complex nature of subglacial sediment transport and the transitions between supply- and transport-limited regimes. Equivalent values of water input and sediment transport capacity below the glacier result in simulated sediment discharges that vary over orders of magnitude (Fig. 10). In turn, using solely the water discharge or sediment transport capacity (e.g., Eq. 8) fails to consider the changes to sediment availability caused by sediment transport, especially when changes to sediment storage can take place over seasonal to decadal timescales.

6 Conclusions

A two-dimensional subglacial sediment transport model, SUGSET_2D, evolves a till layer in response to changing subglacial hydraulic conditions. The model represents sediment transport in supply- and transport-limited regimes, and sediment and water are routed across the bed in response to changing hydraulic conditions in two horizontal dimensions. The till layer is supplied with sediment either from bedrock erosion or by existing sediment, represented by the initial condition. Model cases utilize geometries and hydrological forcings from a synthetic case and Griesgletscher, an alpine glacier in the Swiss Alps.

The interdependence of a large number of parameters and their interactions with one another – for instance, sliding and erosion (Eqs.  12 to 15) – at the very least point to the complexity of sediment transport in the subglacial system. Furthermore, the model's limited representation of the magnitude of interannual variability in the Griesgletscher simulation, from 2011 to 2017, points to processes not completely represented in this application of the model. This misfit could come from poorly constrained parameters and external factors, such as model inputs that may limit the model's accurate representation of sediment discharge observations. These include interannual variability of glacier velocity and thus bedrock erosion, changing glacier topography that routes water to different patches of the glacier bed over time, and routing of water to the glacier bed.

Additional insights into subglacial erosion and sediment transport processes over decadal timescales can be gained from more sophisticated parameterizations of bedrock erosion and subglacial hydrology. Even so, the foundational processes of the model presented here should be considered when examining subglacial sediment transport processes at seasonal to decadal scales. These processes include (1) fluvial transport of subglacial sediment across a glacier's bed in two dimensions in supply- and transport-limited regimes, (2) spatially distributed bedrock erosion or sediment production, and (3) variable water routing in response to changing melt and hydraulic conditions. It is our hope that the model will be applied in the context of field observations to evaluate and isolate subglacial processes controlling sediment discharge from glaciers as they change.