Modeling the spatially distributed nature of subglacial sediment transport and erosion

. 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 transport sediment in one dimension, usually along a glacier’s ﬂow line. Such models have proven useful in describing the formation of landforms, the impact of sediment transport on glacier dynamics, and the interactions between 5 climate, glacier dynamics, and erosion. However, in one dimension, 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. In turn, the geoscience community needs modeling frameworks that describe subglacial sediment discharge in two spatial dimensions ( x and y ) 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 10 glaciers. Experiments performed using an idealized alpine glacier demonstrate the heterogeneity in sediment transport and bedrock erosion below glaciers. The experiments show a non-linear increase in sediment discharge following increased glacier melt. We also apply the model to a real alpine glacier. We compare simulations with annual observations of sediment discharge measured from Griesgletscher in the Swiss Alps. SUGSET_2D reproduces the year-to-year sediment discharge pattern measured at the glacier terminus. The model’s ability to match the data depends greatly on the sediment grain size parameter, which 15 controls subglacial sediment transport capacity. Smaller grain sizes allow sediment transport to occur in regions of the bed with reduced water ﬂow and channel size, effectively increasing sediment connectivity into the main channels. Model outputs from both cases show the importance of considering heterogeneities in water discharge and sediment transport in both the x − and y − dimensions.


Introduction
Increasing glacier ablation perturbs the ways that 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 downstream social and earth systems (Milner et al., 2017;Li et al., 2021).In turn, predictive models are needed to understand the response of these systems to glacier retreat.In alpine environment, increased sediment discharge transported (e.g.Lane et al., 2017;Delaney et al., 2018b).As a result, reducing the problem 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).The one-dimensional models to date 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, describing subglacial sediment transport inherently lends itself to a discretization of bedrock erosion, sediment transport, water flow, and sediment availability in the downstream and transverse dimensions (x and y).
In this manuscript, we present SUGSET_2D, a new two-dimensional subglacial sediment transport model.The model implements subglacial the sediment transport and bedrock erosion processes.We implement a routing scheme that transports sediment in x and y based on the local hydraulic potential gradient.Synthetic cases demonstrate the model's ability to reproduce known processes and yield insight into the spatially-distributed processes responsible for subglacial sediment dynamics.
We also apply the model to a real alpine glacier, Griesgletscher in Switzerland.The model was run with hydrology and topography data from the glacier and measured sediment discharge data are used to validate the model.Through these experiments, we explore the importance of two-dimensional sediment connectivity in the subglacial environment.

Model Description
The model presented here implements a hydraulic model and sediment routing scheme that translates the one-dimensional subglacial sediment transport model presented in Delaney et al. (2019) to two dimensions.In this section, we review the underpinnings of the model presented in Delaney et al. (2019), describe the routing scheme, and outline its numerical implementation in two dimensions.

Hydraulic Model
SUGSET_2D requires a hydraulics model as a means to route sediment and water through the subglacial environment.The hydraulics model is also needed to evaluate the sediment transport capacity of this water, based upon the hydraulic gradient, channel size and water flux (Table 1, Section 2.2; e.g.Walder and Fowler, 1994;Alley et al., 1997).The hydraulic model is based on the premise that subglacial water flows along the hydraulic potential gradient and that the weight of ice pressurizes water at the bed (Shreve, 1972).We simulate characteristics of an Röthlisberger-channel without explicitly describing properties such as creep closure and pressure melt of channel walls.
The hydraulic gradient of a subglacial channel Ψ (at a certain location and time) can be determined with a known hydraulic diameter D h and water discharge Q w .The hydraulic gradient can then be determined using the Darcy-Weissbach equation for fluid flow through a pipe the density of water is ρ w ,the Darcy-Weisbach friction factor is f r , and the channel geometry is accounted for by s (Hooke et al., 1990) s can be represented as 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.The channel width w c is given by where S is the cross-sectional area of the channel given by To approximate the hydraulic diameter D h , we assign a representative water discharge Q * w to Q w , by taking a characteristic water discharge over a certain time period prior (hours to days).We assume that the hydraulic diameter of the channel results from this characteristic water discharge we call the source percentile (c.f.de Fleurian et al., 2018;Delaney et al., 2019;Nanni et al., 2020).The response time and source percentile remain 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, where as prolonged melt would increase the hydraulic diameter.
Water storage is not allowed, such that Q * w and Q w comprise the total amount of melt water produced upglacier.We then evaluate D h , the hydraulic diameter given Ψ * is a representative hydraulic gradient at overburden pressure, evaluated using the Shreve potential gradient 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 Equation 1 to evaluate the instantaneous hydraulic gradient Ψ.To prevent unreasonable water pressures when Q * w rapidly increases and D h is small, the model limits the minimal cross-sectional area S to 0.5 m 2 .

Till-layer model: bedrock erosion and sediment transport
The model simulates the evolution of a subglacial till layer, which we define as transportable sediment below the glacier due to glacier erosion and fluvial sediment transport.Fluvial sediment transport, in supply-and transport-limited regimes, mobilizes and deposits sediment, adding or removing 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 layer.To represent these processes, we implement the Exner Equation (Figure 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.

∂H ∂t
H is till thickness and t is time ( w is the width of a patch of glacier bed, prependicular to the direction of water flow.Q sc is sediment transport capacity, or the amount of sediment that could be transported 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 that enables smooth transition over the range: . ∆σ is a value below which σ substantially deviates from 1, and reduces sediment mobilization. Condition 8a represents the case where bedrock erosion exceeds sediment mobilization, thus sediment transport exists in a transport -limited regime.Condition 8b impedes mobilization or deposition, 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 maximal change in till height over the model run (∼ 10 cm).This term prevents unbounded sediment accumulation, as the model does not include physical processes to limit sediment deposition, such as reduced channel size in response to infill of sediment (Perolo et al., 2018).Condition 8c allows sediment mobilization to transition between transport-and supply-limited regimes, limiting sediment mobilization to sediment production term ṁt (see below), when H is small.
We calculate sediment transport capacity Q sc using the total sediment transport relationship by Engelund and Hansen (1967), 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.We determine the shear stress through the Darcy-Weisbach formulation: where v = Qw S is the water velocity and S is evaluated in Equation 1.Other sediment transport relationships using shear stress could be exchanged by the model operator (e.g.Meyer-Peter and Müller, 1948).We chose Engelund and Hansen (1967)'s formulation due to the representation of both suspended and bedload transport.Source term ṁt is described as, where H max is a till height beyond which no further erosion, ė may occur.
We chose to use an empirical relationship with sliding velocity u b to describe bedrock erosion, k g is an erodability constant and l er is an exponent, which varies from between 0.66 and 3 (Herman et al., 2021).u b is assumed to be relates to basal shear stress (τ b ; Weertman, 1957) given the following relationship 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) where ρ i is the density of ice, h is the glacier thickness and α is the surface slope of the glacier.
Note that because ṁt is a source term, alternative parameterizations of erosion or basal sliding can easily be exchanged.

Spatial and temporal discretization, and parameters
Here, we describe the numerical implementation of the equations presented above, and in particular the routing scheme that enables a two-dimensional representation of subglacial fluvial and till dynamics.

Numerical implementation and parameters
We use a regular grid to discretize the bed.Spatial discretization must be substantially smaller than characteristic length-scale, l, in Equation 8. We then solve Equation 7 for till height H for given initial and boundary conditions in response to till production ṁt and 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 multistep method (Adams-Moulton) that is well-suited for 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 mn, which varies depending on sediment transport conditions and solver tolerance.Longer time steps occur over periods when glacier melt, and thus sediment transport cease (i.e., winter months).Table 3 presents the numerical parameters used.
We impose boundary conditions on the edge cells so no sediment enters the domain.At outlet cells, a flux of sediment leaves the domain, based on sediment transport conditions.Boundary conditions could also be set to represent processes such as hillslope erosion that route material to the subglacial environment (e.g.Andersen et al., 2015).
Evolving Equation 7requires an initial till height, H 0 , 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 checking 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.

Routing algorithm and implementation
Sediment and water are routed down the hydraulic gradient using a multi-cell routing scheme (Quinn et al., 1991), implemented in a similar way as Bovy et al. (2016), but instead on a regular grid.Sediment and water moves from one cell to another using a steepest-descent algorithm, based upon the hydraulic potential.This routing scheme returns a stack, which contains information about the order of cells to perform the calculations.The model evaluates the hydraulic potential at every time step, first, the flotation fraction for a cell at a given time is calculated by f f = φo φ * , where hydraulic potential φ o comes from Here, Ψ j comes from Equation 1, δ is the area of a cell on a regular grid yielding cell length, n r is the number of receivers that the cell has, and w r is the proportion of hydraulic potential fed by the upstream cell.The operation is executed on a cell by cell basis using the routing scheme above, beginning at the glacier lower elevations and moving up glacier.
We distribute the mean value of f f across the glacier and then implement the routing scheme for the hydraulic potential determined from the Shreve potential as The node ordering algorithm is executed every time step in response to diurnal variations in water pressure and thus variable routing of subglacial water in response to changing hydraulic conditions (e.g.Iken and Bindschadler, 1986;Chu et al., 2016).However, to improve stability during periods of rapidly changing sediment transport conditions, we reorder the stack, based upon the hydraulic conditions to the nearest 6 mn.Smaller solving tolerances increase the computational time due to 1) 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 fill closed basins in hydraulic potential to maintain continuous sediment transport through the domain.The model uses an external algorithm, that contains routes flow and fills basins based upon rasterised values of the hydraulic potential.
Using this routing scheme, we are able to evaluate the water discharge in a cell from melt upstream as where ṁw is a melt water source term and n r is the number of receivers of that cell.The sediment discharge Q s into a cell is like-wise computed as Q s is then used to evaluated the ∇ • Q s given Equation 8and subsequently, the change in till height using Equation 7. In both Equation 18and 19, the operations are conducted given the node ordering information in the stack, such that the flux in to a cell depends on the flux through the catchment above it.

Model Application
We use two cases to highlight model viability under increasingly complex situations.First, we apply the model to a synthetic alpine glacier with synthetic hydrology, based on the Subglacial Hydrology Model Inter-comparison 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, and sediment and water discharge at Griesgletscher in the Swiss Alps.We demonstrate the proficiency of the model by comparing sediment transport model output and data (Delaney et al., 2018a).We also identify some drivers of subglacial sediment discharge in the model from these simulations.

Experiment design
We run simulations using an alpine glacier geometry and hydrological forcing following the SHMIP project experiments (de Fleurian et al., 2018)).The domain is 6000 m on one axis and 1080 m on the other.The resulting geometry approximates the Bench Glacier.The U-shaped bed and variable ice thickness mean that variable hydrologic gradients will occur laterally across the glacier and water can be routed across multiple cells.
To represent hydrology, we implement a simple melt model as in SHMIP ( de Fleurian et al., 2018) where M f = 0.01 m K −1 d −1 is a melt factor.T (z s ) is air temperature T at elevation z s defined as where A a and A d are the annual and diurnal amplitudes, respectively, ∆T is a temperature offset, which is adjusted to control the meltwater input and s day are the number of seconds in one day, s year is the number of seconds in a year and dT dz = −0.0075K m −1 is the air temperature lapse rate.In this case, we route water directly to the subglacial system at the location where the melt occurs, ignoring moulins that concentrate meltwater delivery to the bed.
We run the model for 12 years with a steady climate, then we apply a linear temperature increase for 8 years followed by 10 years of steady temperature at the maximal ∆T .The model is initiated with 10 cm of till across the bed.A spin-up over one year of the initial hydrological forcing is applied for 150 a or an annual change in the till layer height is less than 10 −4 mm a −1 , well below the annual erosion rate in most glacierized catchments (Hallet et al., 1996).

Model outputs and findings
Simulations show that over seasonal timescales, sediment discharge increases at the onset of melt and decreases shortly thereafter, prior to the maximum amount of water discharge that occurs each melt season (Figure 4).Daily-averaged sediment discharge decreases until the very end of the melt season when sediment discharge increases slightly again (Figure 6).This occurs when water stops flowing during the night allowing sediment to accumulate in the channels from bedrock erosion.Increased sediment discharge 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 limits transport sediment.
Increases in sediment discharge at the onset of melt have been observed for real glaciers (Willis et al., 1996;Swift et al., 2005;Riihimaki et al., 2005;Delaney et al., 2018b) and reproduced in the one-dimensional version of this model (Delaney et al., 2019).However, in SUGSET_2D, larger diurnal increases in sediment discharge occur near peak daily melt because the area of flowing water expands under the glacier.As a result, increased sediment transport may occur in cells with substantial sediment when hydraulic conditions permit, then abandoned them when water is routed to another part of the glacier bed.This allows sediment to be stored in these cells, until the hydraulic conditions return and increased sediment transport may return.
Such a process is difficult to represent in a one-dimensional model, where the many of the cells could be represented together.
Over the course of the simulation, the mean till height decreases and more sediment is expelled from the glacier as the system adjusts to the change in climate (Figure 5).This occurs despite the spin-up threshold of 10 −4 mm a −1 change in till height per year, highlighting the difficulty in achieving a true equilibrium between bedrock erosion and sediment transport (Figure 4,c).
Till height decrease accelerates following the onset of increased melt at year 12.With the new steady climate reached at year 20, the annual quantities of sediment discharge began to decrease.This occurs as the system approaches a stable relationship between sediment transport and bedrock erosion.
Interestingly, the model recreates "first-flush" events of increased sediment discharge early in the melt season (Figure 6), followed by decreased sediment discharge (Swift et al., 2005;Delaney et al., 2018b).This seasonal evolution in sediment discharge is attributed to increased access to subglacial sediment early in the season, followed by decreased access as flow becomes increasingly channelized (e.g.Willis et al., 1996;Swift et al., 2005).As melt begins in the spring, sediment discharge increases largely due to increased sediment transport high on the glacier (Figure 4).However, maximum values of sediment discharge during "first-flush" events decrease because sediment transport also occurs over that time.By the end of the simulation, winter transport is large enough to prevent the till layer from growing.This maintains a sediment reservoir available for transit when melt increases.Note that the model does not couple ice dynamics to sub-glacial hydrology, so erosive potential, water discharge and the subglacial area will decrease as well in response to the changing climate.Additionally, increased subglacial water discharge could enhance sliding, and thus erosion, may occur following the onset of melt in the spring (Ugelvig et al., 2018).
For the cases described above, bedrock erosion relies only on driving stress and till thickness.Sliding and bedrock erosion did not vary seasonally (Figure 6 a, b, Section 3.1).This causes sediment to accumulate during the winter months, which subsequently provides ample material for transport when melt increases in the spring.The model's bedrock erosion scheme is most applicable to land-terminating glaciers and over long-time scales when driving stress likely exerts a primary control on glacier sliding (e.g.Weertman, 1957;Gimbert et al., 2021).However, by coupling subglacial hydrology to erosion Ugelvig et al. (2018) shows that erosion varies seasonally and abrasion largely occurs solely during the summer months.Additionally, in the case presented above, sediment production occurs primarily near the glacier front, where driving stress, and thus sliding, is highest (Figure 3, a).
To test the effects of spatially variable erosion and the role of hydrology, we present two additional cases to supplement the alpine glacier case above, ORIGINAL.The first case, SEASON, simulates bedrock erosion by increasing sliding during the summer months (e.g.Iken and Bindschadler, 1986), the same erosion relationship is applied as the case as Section 3.1.In this case however, erosion only occurs when the amount of water input substantially exceeds the background basal melt input rate, that is present in the winter.In the second case, CONST, bedrock erosion remains constant over the entirety of the glacier at a rate of 1 mm a −1 .
The ORIGINAL case discharges over 11620 m 3 of sediment per year, while the SEASON case discharged only 60% of that value due to the absence of bedrock erosion during the winter months.The CONST case discharged 7320 m 3 of sediment over the year.CONST's quantity of sediment discharge results in roughly 1 mm a −1 erosion rate due to decreased erosion efficiency with till height (Equation 12 and the limited portion of the bed over-which sediment transport occurs (Figure 5).
Over the three cases, sediment discharge increases at the onset of melt and substantially decreases by the end of the melt season due to sediment exhaustion.In ORIGINAL (Figure 6 a, b), more sediment discharge occurs compared to the alternate cases (SEASON and CONST).The increased sediment discharge in ORIGINAL is due 1) to the prolonged period over which bedrock erosion occurs adding more sediment to the layer and 2) that bedrock erosion occurs low on the glacier where much sediment transport takes place, compared to the CONST case.The peak sediment discharge in CONST (Figure 6 e, f) occurs slightly earlier in the season, due to the increased amounts of sediment on the lower glacier margins.

Experiment design
We also simulate Griesgletscher in the Swiss Alps using topographic data from (Delaney et al., 2019).Hourly water discharge from the glacier was modeled in Delaney et al. (2018a).Here, we use the discharge time series from 2009-2017.Subglacial sediment discharge from the glacier was determined for four different time periods since fall 2011 by differencing bathymetry maps (Delaney et al., 2018a).To estimate surface melt across the glacier with respect to elevation, we use, γ is the mass balance gradient and z 0 s represents the glacier's lowest elevation.ḃ0 represents the melt rate at the glacier's lowest extent.ḃ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 ; a primary control on fluvial transport of subglacial sediment), sliding rate factor (B; a control on bedrock erosion), and the initial till height condition (H 0 ; to approximate the effects of existing quantities of sediment below the glacier).100 simulations were run with randomly selected parameters from a uniform distribution.No spin-up was applied in this case, because of the wide range of H 0 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 was done in the previous cases, the maximum value was applied with an upper limit of 1.
We only considered model outputs resulting in a perfect rank correlation across the four data collection periods and an error less than 80, 000 m −3 .For the case presented below, we show the simulation with the lowest absolute error between model output and the sediment transport data.

Model outputs and findings
The model reproduces the interannual variability in sediment discharge from the Griesgletscher.The absolute error between the model and the measurements is roughly 62, 600 m 3 .The model represents the last three measured yearly sums of sediment discharge well, but it has trouble reconciling the elevated sediment discharge in 2012 and 2013 (Figure 7).This suggests that processes not adequately represented in the model are responsible for the increase in sediment transport, such as activation of new patches of the glacier bed or the relocation of channels (e.g.Zechmann et al., 2020), potentially due to changes to glacier surface topography that cause alternative flow paths below the glacier.Furthermore, glacier sliding, remains constant over the model run, in turn, the results do not explicitly account for seasonal or interannual variability in bedrock erosion (e.g.Herman et al., 2015).
The error from this parameter search is slightly less than half of the 131, 300 m 3 total sediment discharged from the Griesgletscher over this time period (Delaney et al., 2018a), and the error is slightly more than the 58, 300 m 3 from the best model run of the one-dimensional model in Delaney et al. (2019).However, in contrast to the ensemble model runs in Delaney et al. (2019), this model's ability to reproduce the validation data largely depends on the grain size parameter, D m .Compared to Delaney et al. (2019), the sliding parameters and initial condition parameters (B and H 0 ) have a minimal influence here when tuned to the data Figure 7.The dependence on grain size for SUGSET_2D results from the subglacial sediment connectivity parameterized in this two-dimensional version of the model.The channelized nature of flow means that sediment transport may only occur over a relatively narrow patch of the glacier bed (Figure 9).As sediment grain size decreases, sediment from locations of the glacier bed with relatively small water velocity and discharge can more easily be transported to the main glacier channel and be expelled from the glacier.In a one-dimensional model, sediment access occurs over the entire width of the glacier bed.Thus, the bedrock erosion or sediment production term (largely controlled by sliding rate factor B) represents this process, and increased sediment production results in greater connectivity.
The size and shape of the subglacial channels contribute to the discharge of sediment, as well.The sediment transport due to the velocity of subglacial water is limited by the channel width in smaller channels (w c , Equation 10).For this reason, sediment exhaustion occurs mainly in main channels, where channel widths are sufficiently large to allow substantial sediment transport (Figure 3).Conversely, sediment persists in patches of the glacier bed where water velocity could be high, but insufficient channel size effectively reduces sediment transport capacity.Increasing the friction factor f i increases the area of the glacier bed over-which water with substantial velocity flows in SUGSET_2D.Thus the model has trouble capturing interannual variability because sediment exhaustion does not occur over a substantial portion of the glacier bed.
The value of B, from the parameter search, results in a average of 39 m a −1 of glacier sliding across the glacier bed, and the range of values for B in the parameter search result in mean sliding velocities between 14 ma −1 and 70 ma −1 .Yet, due to the low dependence of sediment transport from the glacier on B, other values could perform well but are not captured in the relatively small number of model runs herein.However, because sediment production decreases with till height (Equation 12), sediment production 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, even though the sliding and erosion parameters applied here are within a well constrained range.At the same time, the limited spatial extent of glacier erosion and sediment transport points to a need to evaluate the precise location of bedrock erosion and the impact of subglacial till layers on bedrock erosion in future research.
The best performing model run shows strong temporal variability in sediment discharge (Figure 8), with water discharges from the glacier above roughly 2 m 3 s −1 responsible for much of the sediment transport.Despite the strong dependence on grain size and fluvial transport of sediment in the inversion, sediment transport capacity Q sc remains roughly an order of magnitude higher than sediment discharge (Q sc ).The steep section of the glacier experiences sediment depletion over the model run, as do several channels near the over-deepening and high on the glacier (Figure 9 c d).On some parts of the upper glacier, bedrock erosion in the absence of substantial sediment transport is visible.With changing melt patterns or evolving glacier hydraulic gradients, this sediment could be mobilized and increase sediment discharge down glacier.
The lack of knowledge regarding the spatial distribution of subglacial sediment makes selecting an initial value of H difficult.
The slow rate of basal erosion means that an equilibrium between fluvial sediment transport and bedrock erosion will likely take centuries to attain, if such an equilibrium may even exist in light of variable climatic, and thus glacier, conditions.Should an equilibrium exist, it is probably outside of a feasible computational time given current processing speeds (e.g.Herman et al., 2018;Delaney and Adhikari, 2020).
In addition to selecting an initial value of H, we also limit the thickness at which the till must stop accumulating (Equation 8b, H lim ) due to changes in the hydraulic potential caused by channel infill.We assume that this value is on the order of tens of centimeters (Table 2), based upon available observations (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), on the same order of H lim to improve computational time.Additionally, the model does not consider the interactions between fluvial transport of sediment and debris concentrations in subglacial ice, which may be important for sub-glacial sediment transport (e.g.Ugelvig and Egholm, 2018).SUGSET_2D also contains 20 parameters (Table 2 and 3).These parameters have only 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).
The routing method we use assumes that water solely flows in response to the hydraulic potential (Section 2.3.2).Our parameterization does include the impact of a channel's size on the hydraulic potential.It also 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 (Figure 3; e.g.Werder et al., 2013;Zechmann et al., 2020).Furthermore, a response time of the subglacial channel is chosen prior to simulations to improve computational time, compared to a more sophisticated representation or processes in an R-channel (e.g.Röthlisberger, 1972).

Implications
Results of both the one-dimensional model (SUGSET; Delaney et al., 2019) and SUGSET_2D highlight the importance of simulating the spatial heterogeneities in bedrock erosion, sediment availability, and sediment transport capacity.Yet, in SUGSET, only the till layer (e.g.Equation 8c) and variations in sediment access along the glacier flow line impact sediment transport.In SUGSET_2D , sediment access and transport is 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 their impact on subglacial sediment dynamics (Figures 5, and 9).
Large diurnal and seasonal fluctuations in sediment transport in the synthetic alpine glacier case result from diurnal and seasonal variations in water routing and thus increased sediment availability because sediment transport only occurs over a patch of bed for a short amount of time (Section 3.1).For instance here, diurnal fluctuations in sediment discharge in the middle of the season can be 50% above the mean value, which aligns more closely with some field observation of sediment discharge (e.g.Swift et al., 2005;Delaney et al., 2018b) compared to SUGSET (c.f.Delaney et al., 2019).Furthermore, the results show that the location of bedrock erosion, processes in the till-layer, and the timing of melt all play an important role in the quantity of sediment discharge and the peak sediment discharge that is reached.
In the final case, we compared model runs across a parameter space to sediment discharge data from Griesgletscher in the Swiss Alps (Section 3.2).These results depended solely on sediment grain size compared to the initial till condition or bedrock erosion (Figure 7).Grain size is a strong control in the SUGSET_2D because it modulates how easily sediment patches only accessed by sub-glacial flow during the melt season are mobilized.This process cannot be considered in a one-dimensional model, though it is important even in this relatively small and shallow alpine glacier.These results show that connectivity between subglacial channels and distal sediment patches is a strong control on sediment discharge from the subglacial system.
The connectivity between the main channels and distal sources of sediment could be through the transport of small sediments as applied here, but may also occur through other processes not considered in the model, such as till deformation (e.g.Damsgaard et al., 2020).
Lastly, the model demonstrates the complex nature of subglacial sediment transport and the transitions between supply-and transport-limited regimes.Sediment discharge depends not only on hydrology but also on the sediment availability.Equivalent values of water input and sediment transport capacity below the glacier result in simulated sediment discharge that vary over orders of magnitude (Figure 10, a).In turn, using solely the water discharge or sediment transport capacity (e.g.Equation 10) fails to consider the changes to sediment availability caused by sediment transport, especially when changes to sediment storage can take place over seasons to decades.Finding ways to evaluate these difficult to measure parameters could be key to improving our understanding of subglacial sediment transport.

Conclusions
This manuscript presents a two-dimensional subglacial sediment transport model, SUGSET_2D, that evolves a till-layer in response to subglacial hydrology.Model cases utilize geometries and hydrological forcings from a synthetic and a real alpine glacier.The model captures sediment transport in supply-and transport-limited regimes.Results from both cases point to the need to quantify the spatial distribution of subglacial sediment and water when simulating sediment discharge expelled from glaciers.Model outputs reproduce many observed subglacial sediment processes.
Despite the model's ability to reproduce observations, it relies on a large number of poorly constrained parameters.To our knowledge, only one study has quantified till thickness at a single point below a glacier (Truffer et al., 2000).These observations are limited due to the difficulty of making direct observations at glacier beds.The initial till height in the model must be chosen carefully because the system can remember this initial condition for centuries.The interaction of bedrock erosion and fluvial sediment transport also leads hysteresis in the system.Two-dimensional sediment transport models can represent more observed characteristics of subglacial sediment discharge compared to one-dimensional models.SUGSET_2D routes water and sediment using the Shreve potential and a spatially uniform flotation-fraction that evolves in time in the real glacier case (e.g.Section 3.2).Future work may consider using a coupled model of channelized and distributed drainage networks (Hewitt, 2013;Werder et al., 2013).Increasing the sophistication of the subglacial hydrology model may better evaluate the locations of high potential sediment transport.Such models could even be run offline if the operator assumes, as we do, that rates of change in till height are small compared to the evolution in cross-section of the subglacial conduit.
Our simulations highlight that increased glacier melt does not necessarily result in commensurate changes to sediment discharge unless new previously inaccessible subglacial sediment patches are accessed by meltwater.Additionally, results demonstrate the role of spatially varying water routing and lateral sediment connectivity in subglacial sediment discharge.Further efforts should constrain the role of climate change on glacial dynamics, erosion and sediment transport.Further modeling and observational studies are needed to better constrain the timescales over-which these processes occur in response to climate change.Bedrock Till    .Water discharge, an input modeled for Griesgletscher in (Delaney et al., 2018b), and sediment discharge, output of the model, from Griesgletscher (a).Below is modeled outputs of sediment transport capacity and average till height (b).Note that sediment discharge capacity is roughly one order of magnitude larger than sediment transport discharge.Additionally, the reduction in till height H through this model run shows that sediment is transported from the glacier bed at a greater rate than it is produced.

Figure 1 .
Figure 1.Cartoon of erosional and sediment transport processes considered in model below image of Griesgletscher in 2016.Bedrock erosion scales with sliding speed (us) and adds material to the till layer with thickness H, while water (Qw) transports sediment (Qs) fluvially, if sediment persists in that location of the glacier bed and fluvial transport conditions are sufficient.

Figure 2 .Figure 3 .Figure 4 .Figure 5 .
Figure 2. Illustration of terms in Equation 7, detailing the layers of bedrock, till, water and ice.Characteristics of the subglacial channel are also noted, but shown in one dimenstion for clarity.

Figure 6 .Figure 7 .
Figure 6.Annual response to different erosion patterns across the glacier, although diminishing bedrock erosion with respect to till height is still in place given Equation 12. (a,b) Conventional model setup, where sediment is produced year around.This results in the peak amounts of sediment discharge occur in scenario 1 (a,b) where large amounts of sediment accumulated at the glacier terminus during the winter months.(c,d)Equivalent setup to previous, except sediment is only produced in summer months, when water is present at the glacier bed.Thus, till height remains constant over the winter months.(e,f) Steady erosion of 1 mm a −1 across the entire glacier, with no spatial or temporal variability in sediment production.Yet, the different bedrock erosion scenarios each demonstrate increased sediment discharge at the onset of melt and subsequent exhaustion over the course of the season.
Figure8.Water discharge, an input modeled for Griesgletscher in(Delaney et al., 2018b), and sediment discharge, output of the model,

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

Table
).The first term represents fluvial sediment transport processes, where ∇•Q s represents sediment mobilization in either supply-or transport-limited regimes.The second term captures bedrock erosion processes, where ṁt is a bedrock erosion rate.

Table 2 .
Physical model parameters and constants

Table 3 .
Numerical model parameters