Efficacy of bedrock erosion by subglacial water flow

Bedrock erosion by sediment-bearing subglacial water remains little-studied; however, the process is thought to contribute to bedrock erosion rates in glaciated landscapes and is implicated in the excavation of tunnel valleys and the incision of inner gorges. We adapt physics-based models of fluvial abrasion to the subglacial environment, assembling the first model designed to quantify bedrock erosion caused by transient subglacial water flow. The subglacial drainage model consists of a one-dimensional network of cavities dynamically coupled to one or several Röthlisberger channels (R-channels). The bedrock erosion model is based on the tools and cover effect, whereby particles entrained by the flow impact exposed bedrock. We explore the dependency of glacial meltwater erosion on the structure and magnitude of water input to the system, the ice geometry, and the sediment supply. We find that erosion is not a function of water discharge alone, but also depends on channel size, water pressure, and sediment supply, as in fluvial systems. Modelled glacial meltwater erosion rates are 1 to 2 orders of magnitude lower than the expected rates of total glacial erosion required to produce the sediment supply rates we impose, suggesting that glacial meltwater erosion is negligible at the basin scale. Nevertheless, due to the extreme localization of glacial meltwater erosion (at the base of R-channels), this process can carve bedrock (Nye) channels. In fact, our simulations suggest that the incision of bedrock channels several centimetres deep and a few metres wide can occur in a single year. Modelled incision rates indicate that subglacial water flow can gradually carve a tunnel valley and enhance the relief or even initiate the carving of an inner gorge.


Introduction
Textbook descriptions of glacial erosion detail mechanisms of abrasion and quarrying, but mention erosion by subglacial meltwater as a potential, unquantified, additional incision mechanism (e.g.Bennett and Glasser, 2009;Anderson and Anderson, 2010).This imbalance reflects the deficiency in our understanding of the latter.In fact, subglacial meltwater loaded with sediment has been inferred to carve metre-scale channels in bedrock (e.g.Glasser and Bennett, 2004), often called Nye channels (N-channels; Weertman, 1972) and kilometre-scale tunnel valleys (e.g.Ó Cofaigh, 1996;Glasser and Bennett, 2004;Jørgensen and Sandersen, 2006;Dürst Stucki et al., 2010;Kehew et al., 2012), and more recently its role has been invoked as a necessary mechanism in the carving and deepening of inner gorges (Dürst Stucki et al., 2012;Jansen et al., 2014).

F. Beaud et al.: Erosion by subglacial water flow
In numerical models of glacial erosion, sediment transport by subglacial water flow is usually neglected (e.g.MacGregor et al., 2000MacGregor et al., , 2009;;Brocklehurst and Whipple, 2002;Tomkin and Braun, 2002;Anderson et al., 2006;Herman and Braun, 2008;Egholm et al., 2009Egholm et al., , 2011aEgholm et al., , b, 2012;;Herman et al., 2011;Beaud et al., 2014), often under the assumption that sediment is removed instantaneously.Bedrock erosion resulting from subglacial water is likewise neglected.A better understanding of the processes that lead to the formation of tunnel valleys or inner gorges is also important for the evaluation of deep geological repositories for nuclear waste in regions facing a potential future glaciation (e.g.Iverson and Person, 2012).
Subglacial water flows through two main types of drainage systems: distributed and channelized.In numerical models (see Flowers, 2015, for a review), a distributed drainage system is typically represented by a network of connected cavities (e.g.Lliboutry, 1968;Iken, 1981;Kamb, 1987;Schoof, 2010;Bartholomaus et al., 2011), a macroporous sheet of sediment (e.g.Clarke, 1996;Creyts and Schoof, 2009), or a water film (e.g.Weertman, 1972;Le Brocq et al., 2009).These representations reflect field observations of an increase in water pressure with discharge (e.g.Iken and Bindschadler, 1986;Nienow et al., 1998).A channelized drainage system is most often described by a single Röthlisberger channel or a network thereof (e.g.Röthlisberger, 1972;Nye, 1976;Flowers et al., 2004;Kessler and Anderson, 2004;Schoof, 2010;Hewitt et al., 2012).Water velocities are relatively high in the channelized system and, under steady-state conditions, water pressure decreases with increasing discharge (Röthlisberger, 1972).Conduits carved both in sediment and ice, so-called canals (e.g.Walder and Fowler, 1994;Ng, 2000;Kyrke-Smith and Fowler, 2014), have properties closer to those of a distributed system as pressure often increases with discharge (Walder and Fowler, 1994).In winter, the distributed drainage system evacuates most basal water, and water pressures tend to be relatively high.As surface melt becomes significant, the water input becomes too large for the distributed drainage system alone, water pressure increases, and R-channels start to form.Once an efficient drainage system is established, meltwater is routed relatively quickly downstream and baseline water pressures are generally lower than in winter, with large daily fluctuations.As surface melt decreases, channels are reduced in size by ice creep and may eventually close.We hypothesize that this cycle may have a significant effect on bedrock erosion by subglacial meltwater flow (e.g.Willis et al., 1996;Swift et al., 2005).
Three main processes produce bedrock erosion in rivers: abrasion, macro-abrasion, and quarrying (e.g.Whipple et al., 2000Whipple et al., , 2013)).Abrasion is the result of particles entrained by the flow (saltating or in suspension) colliding with the bedrock and is governed by the tools and cover effect, whereby particles (i.e.tools) entrained by the flow impact exposed bedrock but can also shield it if they are immobile (i.e.cover; e.g.Whipple et al., 2000Whipple et al., , 2013;;Sklar andDiet-rich, 2001, 2004;Turowski et al., 2007;Lamb et al., 2008;Cook et al., 2013;Scheingross et al., 2014).Macro-abrasion and quarrying both result from dislodgement of blocks and require a relatively high joint density in the bedrock.Macroabrasion occurs when blocks are dislodged as a result of the impact of moving particles, while quarrying is the result of dislodgement by pressure gradients caused by water flow (Whipple et al., 2000(Whipple et al., , 2013)).Over highly jointed bedrock, quarrying and macro-abrasion can produce large canyons under extreme flow conditions (e.g.Bretz, 1969;Lamb and Fonstad, 2010;Baynes et al., 2015).In this study, we limit our analysis to abrasion (Whipple et al., 2000(Whipple et al., , 2013;;Cook et al., 2013) and use well-established models to estimate the erosion due to total sediment load (Lamb et al., 2008) and saltating load only (Sklar and Dietrich, 2004).Alley et al. (1997) suggest that the sediment transport capacity of an R-channel is most affected by changes in water discharge and hydraulic potential gradient and further posit that, in most cases, the hydraulic potential gradient increases downstream (due to steep ice-surface slopes close to the terminus), so that the transport capacity should also increase.Dürst Stucki et al. (2012) calculate an erosional potential of subglacial water based on the hydraulic potential gradient under a valley glacier and find that the erosional potential increases toward the terminus and could explain the deepening of inner gorges during a glaciation (e.g.Jansen et al., 2014).Both studies are, however, quite speculative regarding the processes behind subglacial meltwater erosion.Creyts et al. (2013) were the first to couple subglacial water flow in a distributed drainage system, basal refreezing, and sediment transport in a numerical model to explore the evolution of bed slopes adverse to ice flow close to the terminus.They demonstrate strong feedbacks between sediment deposition/entrainment and hydraulic conditions as well as the importance of daily fluctuations in water input on the sediment flux.
In proglacial studies of seasonal sediment yield, hysteresis between sediment and water discharge is often observed (e.g.Willis et al., 1996;Swift et al., 2005;Mao et al., 2014).It is usually attributed to changes in sediment availability, tapping of new areas of the bed by the developing drainage system, or increased mobilization caused by sudden changes in the subglacial hydraulic system (e.g.Willis et al., 1996).An event during which sediment transport peaks before discharge is usually defined as clockwise hysteresis and is interpreted as the manifestation of an unlimited sediment source (e.g.Mao et al., 2014).The opposite is true for anticlockwise hysteresis.In a study of bedload transport by a proglacial stream in the Italian Alps, Mao et al. (2014) identify a transition throughout the melt season from hysteresis dominated by clockwise events to hysteresis dominated by anticlockwise events.The authors infer that this transition is due to the activation of different sediment sources across the drainage basin.
Evidence of the erosional action of subglacial meltwater flow is widespread in formerly glaciated regions and appears in the form of N-channels, tunnel valleys, and inner gorges.Tunnel valleys are large (a few hundred metres to kilometres wide and up to tens of kilometres long) channel-like features found within the limits of former continental ice sheets (e.g.Ó Cofaigh, 1996;Glasser and Bennett, 2004;Jørgensen and Sandersen, 2006;Dürst Stucki et al., 2010;Kehew et al., 2012) or in Antarctica (e.g. Denton and Sugden, 2005;Rose et al., 2014) in substrata varying from loosely consolidated sediment to bedrock.Their formation is attributed to the action of pressurized subglacial meltwater, and three particular mechanisms have been proposed (e.g.Ó Cofaigh, 1996;Glasser and Bennett, 2004;Kehew et al., 2012): (1) sediment creep toward preferential groundwater flow paths, (2) carving by sediment-loaded subglacial water flow, and (3) erosion caused by large subglacial water floods.
Inner gorges are narrow canyons incised at the bottom of an otherwise U-shaped valley (e.g.Montgomery and Korup, 2011), and are found extensively in formerly glaciated mountain ranges like the Alps (e.g.Montgomery and Korup, 2011;Dürst Stucki et al., 2012).The origin of inner gorges was originally entirely attributed to postglacial fluvial erosion, although Montgomery and Korup (2011) conclude that such features persist through repeated glaciations instead of being reset by glacial erosion.For example, Valla et al. (2010) find that fluvial incision rates on the order of a centimetre per year occurred during at least the past 4000 years in a gorge in the French Western Alps.Recently Dürst Stucki et al. (2012) showed that pressurized water flow is necessary to explain the longitudinal profile of an inner gorge in the foothills of the Alps, and Jansen et al. (2014) infer from cosmogenic nuclide exposure that the timing of carving of seven inner gorges in the Baltic Shield matches the timing of glacial cover.Inner gorges are therefore most likely the combined product of fluvial erosion during interglacial periods and subglacial meltwater erosion during glacial periods, although the importance of fluvial vs. glacial conditions is probably dependent on surrounding topography.Nye channels, tunnel valleys, and inner gorges share some characteristics suggesting a common genetic origin, although the specific combination of processes responsible for their evolution may differ slightly.
We implement a one-dimensional (1-D) model of subglacial water flow in which a network of cavities is dynamically coupled to one or a few channels.We then compute the shear stress exerted on the bed and use it to compute transport stages and bedrock erosion rates by abrasion caused by particle impacts.We compare the results of models that treat erosion by saltating particles only (Sklar and Dietrich, 2004) and erosion by both the saltating and suspended particles (Lamb et al., 2008).We use the word "erosion" only to describe the carving of bedrock, while "transport", "mobilization", or "entrainment" refer to the movement of unconsolidated sediments.
We first perform steady-state simulations with the channelized drainage system alone to demonstrate basic model behaviour.We investigate the role of ice geometry, surface melt, and sediment supply.We then introduce the coupled (distributed and channelized) hydraulic system and water forcing with sub-seasonal fluctuations to test the importance of transients in the subglacial drainage system.Finally, we discuss the implications of our results for the formation of N-channels and consequently of tunnel valleys over bedrock and the persistence of inner gorges through repeated glaciations.Our specific research questions are as follows.
(1) What are the major controls on subglacial meltwater erosion?(2) How important is subglacial meltwater erosion compared to overall glacial erosion?(3) Can ordinary seasonal melt processes lead to subglacial bedrock channel incision (and potentially the formation of an incipient tunnel valley or persistence of an inner gorge)?( 4) What are the implications of the water flow regime in channels for hysteresis and sediment transport?

Subglacial water flow
We use a 1-D model of coupled cavity and channelized drainage styled after previous models (e.g.Werder et al., 2013), with the numerical distensibility parameter of Clarke (2003) and the water exchange term of Hewitt and Fowler (2008).The system of equations describing R-channels admits both a wave-like solution and the solution for water flow which introduces numerical stiffness.To circumvent this problem, Clarke (2003) proposed to treat water as a slightly compressible fluid and he introduced a numerical distensibility parameter, γ .The water exchange term is required as the channelized and distributed drainage systems are otherwise defined as two independent systems (Flowers et al., 2004;Hewitt and Fowler, 2008).

Channelized drainage
As in most subglacial drainage models, we assume the channelized system to be a network of semi-circular Röthlisberger channels (Röthlisberger, 1972;Nye, 1976).The channels are assumed to be linear (along the x coordinate) and parallel to one another if more than one is considered.The conservation of water mass is given by (Clarke, 2003) where S is the cross-sectional area of the R-channel, φ ch is the hydraulic potential, Q ch is the water discharge, ρ w and ρ i are the densities of water and ice respectively, v cc is the rate of creep closure of the channel walls, ḃch is a wa-ter source term to the channel, ϕ is the exchange rate between the distributed and channelized drainage systems, W is the width of the catchment drained by the R-channel, γ is a numerical compressibility parameter (Clarke, 2003), x is the coordinate along the flow line, and t is time.Note that we assume changes in bed topography with time are small enough to be neglected.Therefore ∂φ b /∂t = 0 and ∂p ch /∂t = ∂(φ ch − φ b )/∂t = ∂φ ch /∂t, where p ch is the water pressure in the channel and φ b the hydraulic potential at the bed.The evolution of channel cross-sectional area S with time is where is the dissipation of potential energy, is the energy to maintain the water at the pressure melting point, L is the latent heat of fusion, and thus ( − )/(ρ i L) is the rate of channel opening by viscous heat dissipation.We refer the reader to the Supplement for the description of , , v cc , and Q ch .

Distributed drainage
The distributed drainage system is treated as a network of connected cavities (Kamb, 1987) undergoing turbulent flow.Assuming that water can be stored englacially or subglacially, the conservation of water mass can be written (Werder et al., 2013) e v ρ w g where φ ca is the hydraulic potential in the cavity network, e v is the englacial void ratio, q ca is the water flux in the cavity network, g is the gravitational acceleration, v o and v c are the opening and closure rate, respectively, and ḃca is a water source term to the cavity network.The evolution of the average cavity height h ca with time is We refer the reader to the Supplement for the description of v o , v c , and q ca .

Coupled channelized and distributed drainage
To couple the two drainage systems we use a term to describe water exchange between cavities and channels as a function of their respective pressure differences (e.g.Flowers et al., 2004;Hewitt and Fowler, 2008): where k ex is an exchange coefficient (Table 1) and N ch and N ca are the effective pressures in the channelized and cavity systems, respectively (see Supplement).The system formed by Eqs.(1), ( 2), (3), and ( 4) is solved with the MATLAB pdepe solver.Following Clarke (2003) we write the shear stress exerted on the bed by subglacial water flow as where f b is the Darcy-Weisbach friction coefficient for the bed and u is the depth-averaged water flow velocity defined as u = Q ch /S in a channel and u = q ca /h ca in the network of cavities.Details pertaining to the shear stress calculation are described in the Supplement.

Subglacial water flow model simplifications
We implement the model in a 1-D continuum such that the drainage system is assumed to be fully connected.Distributed and channelized drainage systems are assumed to be saturated.The bed is impermeable and undeformable.We also neglect the routing of supra-and englacial water, the effect of particles or refreezing on water flow constriction, and the feedbacks between sliding speed and water pressure.

Boundary and initial conditions
Atmospheric pressure defines the downstream boundary condition for both the network of cavities and the R-channels: φ ch (x = X L , t) = φ b + 1000 Pa.For most simulations, a noflux condition is applied at the upstream boundary such that Q ch (x = 0, t) = 0, q ca (x = 0, t) = 0. Otherwise, for the simulation in which discharge is constant throughout the domain, we apply a Neumann boundary condition at the upstream node of the channel: Q ch (x = 0, t) = 18 m 3 s −1 .We use a 1-year model spin-up for all the transient simulations as we find no significant difference in results using 1-and 2-year spin-ups.
Earth In river reaches where bedrock fracture density is low (i.e.blocks are larger than 1-2 m in size), abrasion is the primary erosional process (Whipple et al., 2000(Whipple et al., , 2013)).Particles that are mobilized by the flow and move by saltation or in suspension can impact the bed, and the energy released upon impact can cause erosion.To erode its bed, a river therefore requires tools (particles available for transport) and an exposed bed.
2.2.1 Saltation erosion model (SEM) (Sklar and Dietrich, 2004) The rate of erosion caused by a saltating load can be expressed as where I r is the rate of particle impact with the bed, V i is the volume of bed material removed upon impact, and F e the fraction of the bed exposed.The impact rate (I r ) is a function of the number of particles and their saltation trajectories.The more particles and the shorter the saltation length, the higher the impact rate.The volume removed upon impact (V i ) is a function of the energy released on impact and therefore of the particle mass, its speed normal to the bed, and the physical properties of the bed and the particle.The fraction exposed (F e ) determines the areal extent of bedrock vulnerable to erosion and represents the cover effect.All of the quantities described here are given per unit width.Values of V i , I r , and F e are computed as where ρ s is the sediment density, D is the particle diameter, w si is the particle velocity upon impact, Y is the Young's modulus of the bedrock, k v is an empirical rock erodibility coefficient, and σ T is the rock tensile strength; where q s is the sediment supply rate and L s is the hop length of particles; where q tc is the sediment transport capacity.Substituting Eqs. ( 8), (9), and (10) into (7) and simplifying yields Note that when q s /q tc > 1 the erosion rate becomes ėsalt = 0 because the bed is completely covered with sediment.Constants and parameters are listed in Table 2, and the formulations of L s and w si are described in the Supplement.We follow Sklar and Dietrich (2004) in using the Fernandez-Luque and van Beek (1976) formulation of transport capacity: where r = ρ s /ρ w − 1 is the buoyant density of sediment and τ * c the critical value of the Shields stress.The Shields stress τ * is computed as 2.2.2Total load erosion model (TLEM) (Lamb et al., 2008) The total load erosion model (TLEM) is based on the original model of Sklar and Dietrich (2004) and extended by Lamb et al. (2008) to account for the impact of particles in suspension close to the bed.The total load erosion rate is calculated as where A 1 < 1 is a coefficient to account for the upward lifting of particles close to the bed by turbulent eddies (here A 1 = 0.36), c b is the near-bed sediment concentration, w i,eff is the effective impact velocity, and q b is the volumetric flux of sediment transported as bedload.The calculation of these quantities is explained in the Supplement.By analogy with the saltation-erosion model, c b is similar to the impact rate I r , w i,eff to the volume removed upon impact V i , and (1 − q b /q tc ) to the fraction of the bed exposed (F e ).Similarly to the SEM, when q b /q tc > 1 the erosion rate becomes ėtot = 0 because the bed is completely covered with sediment.Constants and parameters used for the TLEM are the same as in the SEM and are listed in Table 2.The major differences between the TLEM and the SEM are that (1) the effective impact velocity w i,eff (TLEM) increases monotonically with transport stage as it accounts for the effect of turbulence causing particle collisions with the bed, whereas w si tends toward 0 at large transport stages (and therefore ėsalt approaches 0) as particles trajectories become parallel to the bed, and (2) the fraction of the bed exposed in the TLEM (1 − q b /q tc ) uses the bedload flux (q b ) as contributing to the cover.This means that as the transport stage increases, the fraction of bedrock exposed in the TLEM is higher than in the SEM, provided a significant portion of the sediment is transported in suspension.

Erosion model simplifications
We calculate erosion rates with a constant particle diameter (D = 60 mm) and assume that this particle size is representative of the median grain size in an R-channel.This choice is consistent with Riihimaki et al. (2005), who report a median grain size corresponding to very coarse gravel in the proglacial area of Bench Glacier, Alaska.We assume www.earth-surf-dynam.net/4/125/2016/Earth Surf.Dynam., 4, 125-145, 2016 that the bedrock has a uniform resistance to erosion, erosion occurs uniformly across the width of a channel, abrasion is the main erosion mechanism, and the cover effect is linear.A linear cover effect means that as long as the bed is partially exposed, newly deposited particles are assumed to cover bedrock rather than previously deposited sediment, hence the exponent of 1 on the "cover" fraction q s /q tc (see Turowski et al., 2007).We also neglect sediment transport and the effect of downstream fining.

Modelling strategy and rationale
The model outlined above involves numerous variables and parameters leading to feedbacks.Applications of the model of erosion by abrasion in rivers (e.g.Sklar and Dietrich, 2004, 2006, 2008;Turowski et al., 2007;Lamb et al., 2008;Chatanantavet and Parker, 2009;Nelson and Seminara, 2011;Egholm et al., 2013) show that the primary dependencies are the transport stage, the relative sediment supply, and the hydraulic potential gradient.Numerical modelling studies of subglacial water flow emphasize the importance of the frequency and amplitude of the water input forcing and of the ice and bed geometry (e.g.Flowers, 2008;Creyts and Schoof, 2009;Schoof, 2010;Hewitt, 2013;Werder et al., 2013;Beaud et al., 2014), although few studies have discussed the shear stress on channel walls (e.g.Clarke, 2003).
In an effort to identify the key variables, parameters, and feedbacks, we start with simple experiments and build up the complexity.The simulations are separated in two subsections: (1) steady-state decoupled simulations of a channelized drainage system only (Table 3) and (2) transient simulations of a coupled channelized and distributed drainage system (Table 4).
In the steady-state simulations (Table 3), we first assess the erosion pattern resulting from a constant discharge along the glacier bed.Then we introduce a water forcing that increases with decreasing ice-surface elevation and test the effect of ice geometry and sediment supply.Additional experiments assessing the effect of water flow through a network of cav-ities, water input, sediment size, and channel wall and bed roughness are shown in the Supplement.In the transient simulations (Table 4) we first analyse the role of a synthetic melt season, then use water forcing following realistic melt seasons, and test the role of ice geometry, channel density, and sediment supply.
As per the theory on which the SEM (Sklar and Dietrich, 2004) and TLEM (Lamb et al., 2008) are based, sediment supply, through the tools and cover effects, exerts a major control on erosion rates and patterns by water flow beneath glaciers.Although it is possible to estimate sediment supply to rivers (e.g.Gurnell et al., 1996;Willis et al., 1996;Kirchner et al., 2001;Orwin and Smart, 2004;Sklar and Dietrich, 2004;Riihimaki et al., 2005;Turowski et al., 2009;Mao et al., 2014), no method exists to quantify subglacial sediment supply or transport.In glacierized catchments, measurements are usually made in the proglacial stream, relatively close to the terminus (e.g.Warburton, 1990;Gurnell et al., 1996;Willis et al., 1996;Orwin and Smart, 2004;Riihimaki et al., 2005;Swift et al., 2005;Mao et al., 2014), where water flow is subaerial and potentially influenced by channel dynamics between the glacier terminus and the measurement station (Warburton, 1990;Orwin and Smart, 2004;Mao et al., 2014).In the steady-state simulations we impose a sediment supply that leads to an interesting and diverse range of simulations illustrating most processes and their feedbacks, whereas in the transient simulations we choose a sediment supply that leads to a sediment yield at the last node equivalent to a few millimetres of glacial erosion per year.
We use a flat bed and a parabolic surface for all but two geometries: STP and WDG (Fig. 1).The STP geometry aims at reproducing the steep front of an advancing ice sheet, while the WDG geometry has a wedge shape that resembles the profile of a thinning and retreating ice sheet margin.We assume that water input is distributed uniformly (no moulins) along the bed except for the simulation S_MOULIN.For almost all simulations we fix the sediment supply per unit width; the total sediment supply therefore increases with channel size.This is similar to the assumption of a uniform till distribution across the width of the glacier in which the channelized water flow sources its sediment.

Results
All quantities related to the erosion model are given per unit width of flow.The cross-sectional area of a subglacial channel (S), and thus its width (W ch = 2 √ 2S/π ), changes more rapidly with distance along the flow path than a typical subaerial river.We need to account for these changes when displaying the results, and therefore introduce three quantities: the total transport capacity Q tc = q tc W ch (m 3 s −1 ), the total erosion computed with the TLEM E tot = ėtot W ch (m 2 s −1 ), Table 3. Summary of steady-state simulations.For simulations in which meltwater input is a function of ice-surface elevation z s , we compute f (z s (x)) = ḃss max × 1 − (z s (x) − z s,min )/z s,max , where ḃss max = 8.5 × 10 −7 m s −1 is the maximum meltwater input rate to the channelized drainage system, and z s,min and z s,max are respectively the minimum and maximum ice-surface elevations.Note that ḃss max = 8.5 × 10 −7 m s −1 corresponds to 7.6 cm of ice melt per day assuming ρ i = 910 kg m −3 .The reference sediment supply used for the steady-state simulations is q s,ref = 3.6 × 10 −3 m 2 s −1 (Table 2).

Simulation
Purpose Forcing Difference from reference run Section S_MOULIN R-channel only Table 4. Summary of transient simulations.In this series of simulations the basal sliding speed is u b = 5 m a −1 and the reference sediment supply is q s,ref = 9.1 × 10 −3 m 2 s −1 (Table 2).Water is fed to the network of cavities rather than the channel, hence ḃch = 0 m s −1 , and the function of surface elevation is that of the steady-state simulations (Table 3).Optimized erosion; q s = 0.6 × q tc 4.2.5 and the total erosion computed with the SEM E salt = ėsalt W ch (m 2 s −1 ).

Steady-state decoupled simulations
We examine basic steady-state behaviour of key model variables such as R-channel cross-sectional area (S), transport capacity (q tc ), impact velocity (w si and w i,eff ), transport stage (τ * /τ * c ), and relative sediment supply (q s /q tc ).In the steadystate simulations, the drainage system is composed of a single R-channel.Simulations are terminated once dependent variables (S and φ ch ) reach steady state.In this series of simulations, unless stated otherwise, we impose the sediment supply q s = 3.6 × 10 −3 m 2 s −1 to produce an interesting and diverse range of model behaviour.

R-channel with constant discharge
In the S_MOULIN simulation we use the reference glacier geometry (Fig. 1, REF) and the water input is imposed at the uppermost node only, i.e. as if a moulin were feeding the system.This permits us to drive the system with a constant water discharge and to analyse the resulting relation between discharge, channel cross-sectional area, velocity, instantaneous erosion rates, and transport capacity patterns.
Although discharge is constant (Fig. 2a), the crosssectional area of the R-channel changes along the profile (Fig. 2a) due to the ice geometry (Fig. 1).Over the first 46 km, the cross-sectional area decreases in response to the steepening hydraulic potential gradient (Fig. 2b), the latter being a function of ice-surface slope (Fig. 1).Close to the terminus (last 4 km of the profile), the ice thins significantly, the hydraulic potential gradient shallows (by a factor of 3), and the cross-sectional area increases.The average water velocity assumes the opposite pattern (Fig. 2b).Because the grain size is kept constant, q tc ∝ (τ * − τ * c ) 3/2 , the rate of sediment transport (Fig. 2c) is amplified relative to the velocity (q tc drops by a factor of 6, while u is reduced by about 50 %); both have maxima at 46 km and decrease sharply near the  3  and 4) and (b) corresponding ice-surface slopes.The geometries REF, IS1300, and IS700 are parabolas such that the surface elevation is given by z s (x) = z s,max √ x, where z s,max is the maximum ice-surface elevation and x ∈ [1, 0].In order to obtain the steeper front in STP we use a cubic root instead of the square root: z s (x) = z s,max x 1/3 .The wedge-like geometry WDG is defined by a straight line with the same elevation change as REF.
terminus.In this simulation the sediment transport capacity is always larger than the supply rate (Fig. 2c), exposing most of the bed (F e > 0.5).
In both the TLEM and SEM, near-bed sediment concentration and impact rate are described as a function of the sediment supply and the hop trajectory of a particle, so they are similar (Fig. 2d).As the velocity increases (km 0-46) the sediment is transported faster and further from the bed and the near-bed sediment concentration and impact rate decrease.
Impact velocities (Fig. 2e) vary depending on the model.In the SEM the impact velocity tends toward zero as the hop length increases and the particles approach transport in suspension, which leads to a local minimum around km 46.The TLEM accounts for the effect of turbulent eddies on the trajectory of particles close to the bed, and thus the effective impact velocity is commensurate with the velocity and peaks around km 46.
Erosion rates (Fig. 2f) in both the SEM and TLEM show a minimum at 46 km and a maximum close to the terminus.The minimum and maximum correspond respectively to the minimum and sharp rise in near-bed sediment concentration or impact rate.Even under a constant discharge, ice-surface slope and channel size produce a peak in velocity just upstream of the terminus.At this peak, the flow has the power to lift particles far enough from the bed to reduce erosion in both the SEM and TLEM.Erosion rates are higher with the TLEM than the SEM due to the difference in impact velocities when the sediment transport regime approaches suspension.For the transport stages we obtain, impact velocities with the TLEM are consistently higher than with the SEM.Moreover, for relatively large transport stages, more of the bed is exposed in the TLEM: when some fraction of the total load travels in suspension, q s > q b .

Ice geometry
The surface slope of a glacier is a first-order control on subglacial water flow (Eqs.1-2).We experiment with thicker ice (IS1300, Fig. 1), thinner ice (IS700, Fig. 1), and steeper icesurface slopes close to the terminus (STP, Fig. 1), all compared to the reference ice geometry (REF, Fig. 1).Shallower surface slopes close to the terminus are tested with a wedgelike geometry of constant surface slope (WDG, Fig. 1).All these simulations (Fig. 3; see Table 3) employ the surfacemelt profile of S_REF and therefore yield nearly identical discharge profiles (Fig. 3a).The water input is a function of the ice-surface elevation (Table 3); hence, the discharge is largest close to the terminus.
The thinner the ice close to the terminus the larger the channel (Fig. 3b) and the lower the transport stage (Fig. 3c).In the case of S_STP, the combination of a particularly steep hydraulic potential gradient and relatively thick ice near the terminus inhibits channel enlargement over the last 5-10 km as observed in other simulations (Fig. 3b).Since the discharge profile is the same for the simulations presented, the smaller the channel, the faster the flow and hence the larger the transport stage.The maximum transport stage for the simulation with the steepest terminus (S_STP, τ * /τ * c ≈ 18.5, Fig. 3c) is almost 4 times larger than that for the simulation with the shallowest terminus (S_WDG, τ * /τ * c ≈ 5.5, Fig. 3c; see Fig. 1 and Table 3).Over the first 25 km of the profile the simulation with a constant ice-surface slope (S_WDG) shows the steepest hydraulic potential gradient, creating comparatively higher transport stages than other models.
Erosion begins around km 4 for S_WDG (Fig. 3d), and around km 19-20 for S_STP and S_700 (Fig. 3d), because for the prescribed sediment supply, the bed becomes exposed for transport stages τ * /τ * c 2.5.All simulations but S_WDG have a local maximum in total erosion (E tot , Fig. 3d) between km 30 and 37 and a local minimum between km 45 and 49.The sediment supply per unit width is constant; therefore, the relative sediment supply decreases as the transport stage increases and thus the erosion rate per unit width decreases when q s /q tc < 0.5 because the number of tools decreases.For simulation S_WDG the increase in channel size compensates for the small decrease in erosion rate per unit width (not shown) because the transport stage remains relatively low (τ * /τ * c < 6).Over the last few kilometres of the profile, total erosion increases again for the simulations in which transport stage drops to moderate values (S_REF and S_1300).Total erosion drops sharply at the terminus for S_700 as the transport stage drops below 2.5.When the transport stage remains relatively high (τ * /τ * c > 15) the number of tools remains low, as does total erosion (S_STP, Fig. 3c and 3d geometry exerts a primary influence on transport stage and erosion patterns via its control on hydraulic potential gradients and channel size.Simulation S_WDG yields the most erosion, despite relatively low transport stages, illustrating the importance of the tools effect (e.g.Sklar and Dietrich, 2006).

Sediment supply
Tools and cover compete so that both a lack and an overabundance of tools hinder erosion.Sklar and Dietrich (2004) and Lamb et al. (2008) have shown that erosion peaks for a given flow regime at an optimum relative sediment supply.We thus investigate how varying sediment supply (q s = 1.8 × 10 −4 to 8.9×10 −2 m 2 s −1 ) affects the rates and patterns of subglacial meltwater erosion, while the subglacial hydraulic regime remains that of S_REF (Fig. 3; Table 3).For sediment supply rates q s ≤ q s,ref × 2 (Fig. 4a) the patterns of total erosion with the TLEM and SEM (Fig. 4b and  4c) show the same features as in Fig. 3: no erosion in the uppermost part of the profile followed by a local maximum around mid-profile, a local minimum around km 46-47 and a sharp rise in the last 3 km.For larger sediment supply q s ≥ q s,ref × 5, the patterns change and have a single maximum only.If the relative sediment supply satisfies q s /q tc > 0.3 (Fig. 4a) the number of tools remains high enough for total erosion to increase with transport stage (τ * /τ * c ) as the increase in channel size compensates for the small drop in erosion per unit width ėtot when q s /q tc < 0.5.
Erosion occurs over 48 km of the bed (S_SSP, q s,ref /20) in both the SEM and TLEM for the lowest sediment supply, whereas the bed is almost completely shielded for the second largest sediment supply (S_SSP, q s,ref × 20) and erosion occurs only between km 46 and 48 in the TLEM (Fig. 4).The SEM treats the whole sediment supply as participat-   3) on total erosion with the TLEM and SEM for the same hydrology as the reference simulation (S_REF, Table 3).The legend applies to all panels.(a) Relative sediment supply (q s /q tc ), (b) total erosion rate computed with the TLEM (E tot ), and (c) total erosion rate computed with the SEM (E salt ).Note the logarithmic scale on the y axis of (b) and (c).
ing to the cover effect, and the bed is shielded as soon as q s /q tc ≤ 1.The TLEM discriminates between transport as bedload and in suspension; therefore, if q b < q s the bed can remain partially exposed even for a relative sediment supply q s /q tc ≥ 1.The erosion window is therefore larger for the TLEM (Fig. 6 in Lamb et al., 2008).The fact that erosion rates computed with the TLEM are higher than those computed with the SEM is inherent to the model formulation (Lamb et al., 2008).An interesting conclusion arising from Fig. 4 is that total erosion is significant over most of the bed (Fig. 4b and  4c) when the sediment supply rate is relatively low (q s ≤ q s,ref × 2; Fig. 4a).Total erosion becomes more localized at higher relative sediment supply rates.Changing the particle diameter (D), instead of sediment supply (q s ), leads to similar changes in relative sediment supply and therefore in erosion patterns (see Supplement).The transport stages calculated are large enough to produce significant differences between the SEM and TLEM for large relative sediment supplies.Hereafter, we focus on the TLEM only in the results as it is more appropriate for the flow conditions encountered beneath glaciers.

Transient coupled simulations
Arguably, most erosion and sediment transport in fluvial systems occur during flood events (e.g.Whipple et al., 2000Whipple et al., , 2013;;Kirchner et al., 2001;Sklar and Dietrich, 2004;Lamb et al., 2008;Turowski et al., 2009;Lamb and Fonstad, 2010;Cook et al., 2013).In glacial environments large daily variations in meltwater input to the glacier bed can be likened to periodic flooding (e.g.Willis et al., 1996).We perform a series of transient simulations with a cavity network coupled to an R-channel (Eqs. 1, 2, 3 and 4) to explore how the transience in subglacial water flow is affected by changes in ice thickness, surface melt, and sediment supply, and how this transience impacts instantaneous and annually integrated erosion rates.In this series of simulations, we choose the sediment supply (q s,ref = 9.1 × 10 −3 m 2 s −1 ) such that the modelled sediment yield ( min(q s (X L )W ch (X L ), Q tc (X L ))dt) corresponds to an inferred basin-wide erosion rate of a few millimetres of erosion per year (e.g.Gurnell et al., 1996;Hallet et al., 1996;Koppes and Montgomery, 2009).

Reference model
The reference model uses a synthetic forcing in the form of water supply to the distributed system ( ḃca , Eq. ( 3); sinusoid with a period of 120 days on which we superimpose daily fluctuations) and the reference glacier geometry (Fig. 1, REF; Table 4, T_REF).The sediment supply rate per unit width is assumed to be uniform along the bed.With this simple test we explore how the transience in water input and response of the subglacial drainage system affect transport stage and erosion rate.Modelling subglacial water flow through coupled distributed and channelized drainage systems has been described extensively in recent literature (see Table 3 in Flowers, 2015), so we omit a discussion of the drainage system itself and focus instead on how subglacial hydrology affects transport stage and erosion.
Since the sediment supply is fixed, transport stage (Fig. 5a) and relative sediment supply (Fig. 5b) are anti-correlated.The time transgression in transport stage and relative sediment supply is a result of the up-glacier incision of Rchannels.Once a channel is well developed, water pressure decreases in the channel, as does water velocity and transport stage (Fig. 5a).Daily fluctuations are only detectable close to the terminus, where a channel is relatively well developed.
Similar to what we find in the steady-state simulations, the largest transport stages are found a few kilometres up-glacier from the terminus.The bed remains shielded over the first 25 km from the ice divide.The time window during which erosion occurs (q s /q tc 1 − 1.5) is longest at km 45 and decreases up-and down-glacier.The size of the R-channel at the terminus (km 50) after day 80 becomes large enough that the transport stage (Fig. 5a) is insufficient to maintain an exposed bed.tential gradient shallows and the discharge decreases so that the bed is exposed for a shorter time.
Erosion per unit width ( ėtot , Fig. 5c) peaks when the relative sediment supply satisfies 0.25 ≤ q s /q tc ≤ 0.4, i.e. at the onset of R-channel formation, even before the peak in transport stage for the lowermost 15 km (km 35-50 in Fig. 5).The peak in erosion is followed by relatively constant values and eventually an abrupt drop.At the terminus, erosion ceases due to low shear stress in the relatively large channel, while upstream erosion ceases due to declining water supply.Note that erosion can occur after the melt season ends (day 120) at km 45 because water remains stored englacially and subglacially (Eq.3).Further up-glacier (km 25-30) the maximum erosion per unit width coincides with the minimum in relative sediment supply, as the latter remains larger than 0.5 and enough tools are available close to the bed.
Patterns of total erosion (E tot , Fig. 5c, thick lines) all peak between day 65 and 75, except at the terminus (km 50).The initial peak in erosion per unit width occurs while the channel is small, and thus total erosion is largely controlled by channel size over most of the record.Simulations with transient meltwater input highlight the role of channel size in controlling transport stages and erosion close to the terminus.

Surface melt
We vary the amount of water reaching the bed to explore the differences between using synthetic and realistic melt records.The realistic melt records come from the ablation area of an unnamed glacier in the Saint Elias Mountains, Yukon, Canada, in 2007and 2008(Wheler et al., 2014).We scale the 2007 melt record (T_2007, Table 4) so that the total volume of water is identical to the synthetic input.The melt time series from 2008 (T_2008, Table 4) is then scaled such that the ratio of 2007 to 2008 melt volumes is preserved.This test is intended to highlight the importance of total melt volume and the temporal structure of meltwater input.
When we apply the realistic forcing from 2007 (T_2007), the transport stage exhibits four to five peaks (Fig. 6a) at km 35, 42.5, and 50.At these three locations, the first peak occurs once enough water is supplied to the bed to form a channel (Fig. 6a; after 45 days at 42.5 and 50 km and after 75 days at 35 km).The subsequent peaks in transport stage (after day 60 at 42.5 and 50 km and after day 80 at 35 km) follow periods of high melt.After day 45, the transport stage Given the prescribed sediment supply rate, the bed is only exposed at transport stages larger than 3.5.Erosion (Fig. 6b; ėtot and E tot ) at km 27.5 and 50 therefore only occurs during peaks in transport stage; at km 20 the bed is always covered.Erosion rate per unit width ( ėtot , Fig. 6b) plateaus at moderate transport stages, thus it remains relatively constant once the bed is partially exposed.On the other hand, total erosion (E tot , Fig. 6b) peaks with transport stage (Fig. 6a).Similar results are obtained for T_2008 (Fig. 6c and 6d), where a different melt time series is employed.While the amplitudes of the fluctuations of meltwater input are a few times larger in T_2008 than T_2007 (Fig. 6a and 6c), the total melt in T_2008 is about 80 % that of T_2007.These realistic meltwater forcings produce episodic variations in transport stage and erosion rate that suggest multi-day fluctuations in meltwater input are important.
These multi-day variations in water input also lead to a succession of channel enlargement events (Fig. 6a and 6c), represented by multiple peaks in transport stage (τ * /τ * c ).On the timescale of several days, creep closure near the terminus is low enough that a channel is sustained between the melt events, leading to an up-glacier migration of relatively large transport stages and integrated erosion.Thus, the pressure in the channel close to the terminus, and hence transport stage, is low.This results in the integrated total erosion ( E tot dt) being about 3 times lower for realistic inputs (Fig. 7, T_2007 and T_2008) than for the synthetic one (Fig. 7, T_REF) at the terminus.For the same total water input (T_REF and T_2007), the realistic melt season produces more erosion averaged over the glacier bed than the synthetic input (8 × 10 −2 mm for T_2007 vs. ∼ 7 × 10 −2 mm for T_REF).

Ice geometry
Studies of sediment yield from glacierized catchments (e.g.Hallet et al., 1996;Koppes andHallet, 2002, 2006;Koppes and Montgomery, 2009) conclude that glaciers are more erosive during retreat than during advance due to the amount of meltwater production.Glacier thinning (or thickening) during a phase of retreat (or advance) will also impact the development of the subglacial drainage system and hence its ability to flush sediments and erode the bed.In this series of model tests we hold the sediment supply fixed and vary the glacier geometry by changing the maximum ice thickness (Fig. 1, T_REF, T_1300, T_700, Table 4), while the water input remains the same.As we have already described the principal mechanisms responsible for fluctuations in erosion rates in previous sections, we now focus on annually integrated erosion.For all ice geometries tested in Fig. 8 (T_1300, T_700 and T_REF, Table 4), significant erosion only occurs downglacier of km 20.The thicker the ice, the further up-glacier significant erosion (both ėtot dt and E tot dt) occurs (up to km 21 for T_1300 and km 29 for T_700).In these tests, thicker ice also means steeper surface slopes (Fig. 1).Since water input is identical for these simulations, steeper surface slopes lead to faster water flow and the possibility of initiating sediment motion further up glacier.At the terminus almost 4 times as much erosion occurs for T_1300 than T_700 because thick ice prevents the growth of a large channel.

Subglacial drainage catchment width of a channel
Hydraulic properties of the distributed drainage system determine the density of channels that form (e.g.Werder et al., 2013).The smaller the channel spacing, the lower the discharge in a single channel.A smaller channel, at equilibrium, yields larger water pressures, so we expect that more water would be evacuated through the cavity network.In this test we fix the total glacier width at 1000 m and allow two, three, or four channels to form such that channel catchment widths (W ) are, respectively, 500, 333, and 250 m (T_W500, T_W333, T_W250; Table 4).If we consider the erosion in a single R-channel per simulation (Fig. 9a), the smaller the drainage catchment width, the smaller the discharge, and the smaller the time-integrated erosion (Fig. 9a; ėtot dt and E tot dt).The feedback causing erosion rate per unit width to decrease at large transport stages (see Fig. 2) is such that even for the simulations where the drainage catchment width is relatively small, erosion rates are comparable ( ėtot dt, Fig. 9a) despite the lower transport stages.The differences are, however, relatively large for the annually integrated total erosion (maximum E tot dt for T_REF is more than twice that of T_W250, Fig. 9a) because of the effect of channel size.
The hierarchy in total integrated erosion is inverted when all R-channels within a fixed glacier width are accounted for (Fig. 9b, between km ∼ 37 and ∼ 48).Once the number of R-channels present is taken into account, integrated total erosion is largest for the smallest channel catchment (almost twice as large for T_W250 than T_REF, Fig. 9b).In this case, numerous small channels therefore produce more erosion than few large ones.
At the terminus (km 50), the simulations with a catchment width per channel smaller than that in T_REF show values of annually integrated erosion ( E tot dt) of about half that of the reference simulation (Fig. 9b).The relatively smaller Rchannels in these simulations remain more pressurized and thus drain less water from the cavity network (Eq.1); the relative discharge in the cavity network near the terminus is therefore larger than in the reference simulation (T_REF), further diminishing transport stage near the terminus (see Figs. 2-3). .Time-integrated erosion per unit width ( ėtot dt) and time-integrated total erosion ( E tot dt) as a function of sediment supply rate (q s ) (Table 4).The first 20 km of the profile are not shown because the bed is alluviated and erosion rates are negligible.

Sediment supply
In the present model the values and patterns of sediment supply are amongst the key unknowns.Most till is produced subglacially (e.g.Sanders et al., 2013) and the amount and size distribution of till depends on the history and patterns of production (quarrying) and comminution (abrasion).We test the sensitivity of erosion rates and patterns to different values of input sediment supply.In two simulations (T_SSP/2 and T_SSP/4; Table 4) the sediment supply rate per unit width is constant in space and time and is taken as a fraction of the reference supply rate in T_REF (Table 4).The largest erosion rate in the SEM occurs when the relative sediment supply is q s /q tc = 0.5.For the TLEM and transport stages τ * /τ * c < 100, the maximum erosion rate is obtained for a relative sediment supply of 0.5 ≤ q s /q tc < 0.8 (Lamb et al., 2008).We determine a ratio q s /q tc close to optimum and examine the resulting erosion rates and patterns (T_SSPOPT; Table 4).This provides us an upper bound on subglacial meltwater erosion rates.The hydraulic conditions in this suite of simulations are that of T_REF (Fig. 5).
Decreasing the sediment supply leads to a decrease in the maximum integrated erosion (Fig. 10, T_SSP/2 and T_SSP/4; ėtot dt and E tot dt) and to the bed being eroded further up-glacier.For a relatively low sediment supply (T_SSP/4), the peak in annually integrated erosion per unit width ( ėtot dt, Fig. 10) is hardly discernible and the peak in annually integrated total erosion ( E tot dt, Fig. 10) is controlled by channel size.In order to estimate the maximum erosion that can occur under the given hydraulic conditions and sediment size, we optimize the sediment supply rate by expressing it as a function of the transport capacity (q s /q tc ≈ 0.6).The resulting patterns of annually integrated erosion ( ėtot dt and E tot dt, Fig. 10, T_SSPOPT) mimic the transport stage patterns (see τ * /τ * c , Fig, 5a), and peak at nearly twice the values of T_REF.

Significance of model simplifications
We have detailed the simplifications and underlying assumptions of the model while describing the model and the strategy; we therefore focus on the potential implications of the most important simplifications.At the onset of the melt season, sliding is expected to accelerate as a response to increased water supply to a distributed drainage system (e.g.Iken, 1981;Hooke et al., 1989;Mair et al., 2003;Anderson et al., 2004;Sole et al., 2011;Meierbachtol et al., 2013;Hewitt, 2013;Hoffman and Price, 2014) which would promote cavity enlargement and water flow through the distributed rather than the incipient channelized drainage system.This sliding feedback alone could produce a small decrease in water pressure (e.g.Hoffman and Price, 2014) and hence a decrease in transport stage.
In this study, we treat only the case of bedrock erosion by abrasion and we neglect the effect of quarrying.Although the latter can lead to erosion rates up to an order of magnitude larger than abrasion, it requires that the bedrock be highly jointed (Whipple et al., 2000(Whipple et al., , 2013)).Quarrying is a two-step process: (1) loosening of blocks around pre-existing cracks (or possibly opening of new cracks) and (2) mobilization and transport of loose blocks (Whipple et al., 2000(Whipple et al., , 2013;;Chatanantavet and Parker, 2009;Dubinski and Wohl, 2013;Lamb et al., 2015).The depth of loose cracks could be related to sediment availability (Chatanantavet and Parker, 2009) and mobilization and transport of quarried blocks scale with the transport stage (Dubinski and Wohl, 2013;Lamb et al., 2015).Therefore we expect that the patterns of quarrying would be similar to the transport stage, yet limited by the thickness of the loosened layer.
We compute erosion with only a single particle size that is assumed to be the median of size of the sediment mixture (e.g.Sklar and Dietrich, 2004, 2006, 2008;Turowski et al., 2007;Lamb et al., 2008;Chatanantavet and Parker, 2009;Nelson and Seminara, 2011).The SEM (Sklar and Dietrich, 2004) was generalized for a grain size distribution by Egholm et al. (2013), a study in which they, however, omit a discussion of the implications of the generalization of the SEM.As for the TLEM, Lamb et al. (2008) suggest that a generalization to grain size distribution would require re-evaluation of some of the equations to account for the interactions between particles of different sizes within the bedload layer.A decrease in median sediment size would probably result in an erosion profile more spread out along the bed and an increase in median sediment size would result in a localization of erosion (see Supplement).The changes in erosion would be quantitatively similar to a decrease in sediment supply (q s ) and thus a decrease in relative sediment supply (q s /q tc ), which strongly controls erosion patterns.
We make the assumption of a supply-limited glacier bed and hence neglect the effect of sediment transport and the in-teractions between sediment thickness and water flow.The mobilization and particularly deposition of sediment affect the flow regime by enlarging or reducing the cross section of flow (Creyts et al., 2013).On a timescale of days, when sediment is mobilized, the cross section of flow is enlarged and could result in a drop in channel water pressure and a corresponding loss of flow strength.The opposite effect, leading to flow strengthening, could occur when sediment is deposited.We do not treat the case of transport-limited conditions, where the channelized drainage system would more closely resemble canals (Walder and Fowler, 1994;Ng, 2000;Kyrke-Smith and Fowler, 2014).To implement sediment transport adequately, it is necessary to improve existing models for subglacial water flow through canals (Walder and Fowler, 1994;Ng, 2000;Kyrke-Smith and Fowler, 2014) with timeevolving effective pressure.Alley et al. (1997), however, argue that in the case of subglacial water flow through canals, the water pressure remains relatively high and very little water would be drained from the distributed system, limiting the capacity of canals to transport sediment.In contrast, the steady-state water pressure in an R-channel decreases with increasing discharge, favouring water flow in the channelized system and enhancing transport and erosion.
Accounting for the production of sediment and the evolution of particle diameter at the glacier bed would also largely influence sediment supply patterns.For example, we can speculate that if the sediment sources are localized in areas of more easily eroded bedrock (e.g.Dühnforth et al., 2010), tools would only be present downstream from these areas.If, instead of fixing the sediment supply per unit width, we fix the total sediment supply (simulation not shown), tools are less available at peak flows, reducing erosion, whereas the cover effect is enhanced for a relatively small channel.We also tested a simple power-law downstream fining function (e.g.Sklar and Dietrich, 2006; simulation not shown).The results were very similar to those obtained with a decrease in relative sediment supply, because particles were smaller than the reference size of 60 mm in the region of the bed where channels form.Another means of obtaining insight into sediment supply rates and patterns would be through the use of a comprehensive model of glacial erosion, i.e. a model encompassing transient subglacial hydrology (between distributed and channelized systems), ice dynamics, glacial abrasion, and quarrying.Such a model is, however, yet to be developed as patterns of glacial erosion remain poorly understood (see Beaud et al., 2014).Finally, subglacial water flow evacuates a significant volume of sediment despite the small area over which R-channels operate and the tendency of these channels to remain stably positioned in association with moulins (e.g.Gulley et al., 2012).The mechanism by which large volumes of sediment are delivered to the channels remains elusive.More work is therefore required to quantify subglacial sediment production patterns.

What are the major controls on subglacial meltwater erosion?
We rank the transient simulations by glacier-area-averaged erosion rate in Fig. 11a.Because we prescribe water input rates sufficient to form a channelized drainage system, it stands out from the model formulation that sediment supply is the most important parameter.A lack or overabundance of tools inhibits erosion.In our results this is shown by the fact that T_SSPOPT (sediment supply optimized for erosion; Table 4) produces the most erosion and T_SSP/4 (smallest sediment supply; Table 4) the least.Changing the ice geometry also leads to a relatively large range of averaged erosion rates as T_1300 (thick ice; Table 4) yields twice as much erosion as T_700 (thin ice; Table 4; Fig. 11a).Larger hydraulic potential gradients in T_1300 cause the shear stress to be large enough over larger portions of the bed to create erosion (Fig. 9).Subglacial drainage catchment width, within the range tested, plays a lesser role than sediment supply or ice geometry although the averaged erosion rate in T_W250 (four channels; Table 4) is ∼ 30 % more than that of T_REF.The fact that T_2007 (realistic melt season from 2007 record; Table 4) produces more averaged erosion than T_REF suggests that an increase in the multi-day variability of the water input enhances erosion.The relations are different for apparent erosion rate (Fig. 11b), here defined as the equivalent thickness of bed material evacuated by the integrated sediment flux ( min(q s W ch , Q tc )dt) at the terminus (km 50).The apparent erosion rate corresponds to the quantity estimated by studies of sediment yield in proglacial channels, lake, or fjords.Relatively large hydraulic potential gradients and relatively thick ice close to the terminus, both of which inhibit R-channel growth, compete against the loss of transport capacity.Therefore the largest apparent erosion occurs for the thickest ice (T_1300, Fig. 3).Interestingly, the lowest drainage density (T_REF) yields more apparent erosion than the highest (T_W250).Discharge through the cavity network close to the terminus increases with R-channel density; the smaller the channel, the larger the water pressure and the lower the pressure gradient between the two systems.This feedback, in addition to the discharge in the R-channel being smaller due to the R-channel drainage catchment size, reduces the transport stage close to the terminus.
The results in Fig. 11b suggest that, despite the increase in apparent erosion that accompanies an increase in meltwater input (e.g. the total melt in T_2007 is about 1.25 times that of T_2008), the thinning associated with the retreat of an ice mass would have a competing effect by decreasing the hydraulic potential gradient (see Fig. 3).The flushing power of subglacial water flow is conducive to the removal of subglacial sediment enabling glacial abrasion and quarrying to be efficient.Our results suggest that the subglacial drainage conditions most favourable for glacial erosion occur where significant surface melt and relatively steep surface slopes occur simultaneously, i.e. during an ice sheet maximum advance or during early phases of retreat.This corroborates the hypothesis of Jørgensen and Sandersen (2006) that some Danish tunnel valleys were excavated during the stagnation of the Scandinavian ice sheet.However, these findings challenge the hypothesis that glaciers deliver more sediment to proglacial areas during retreat than during advance (e.g.Hallet et al., 1996;Koppes andHallet, 2002, 2006;Koppes and Montgomery, 2009), yet more work is required to explore this hypothesis.The lack of flow strength in the upper reaches of the glacier (upstream from km 20 for most simulations) suggests that subglacial sediment in the accumulation area is transported almost solely by entrainment due to sliding at the ice-bed interface.
In the steady-state simulations we find that significant erosion can occur in a network of cavities (see Supplement).In the transient simulations, however, the coupling with Rchannels prevents large shear stresses from developing in the distributed drainage system, and the threshold for sediment motion is not even reached for particles of 1 mm diameter.We thus argue that bedrock erosion in the distributed drainage system is limited unless specific conditions are satisfied, for example a subglacial flood or a surge.

How important is subglacial meltwater erosion compared to overall glacial erosion?
In most literature on modelling landscape evolution by glacial erosion it is assumed that subglacial meltwater efficiently removes sediment from the glacier bed, while its effect on bedrock erosion is neglected (e.g.MacGregor et al., 2000MacGregor et al., , 2009;;Brocklehurst and Whipple, 2002;Tomkin and Braun, 2002;Anderson et al., 2006;Herman and Braun, 2008;Egholm et al., 2009Egholm et al., , 2011a, b;, b;Herman et al., 2011).On the other hand, in formerly glaciated landscapes, erosional features like tunnel valleys (e.g.Glasser and Bennett, 2004;Denton and Sugden, 2005;Dürst Stucki et al., 2010, 2012;Kehew et al., 2012;Jansen et al., 2014) indicate that subglacial water can produce significant bedrock erosion.The results we obtain with our simple ice geometries and water input forcings indicate that the areally averaged bedrock erosion produced by subglacial water flow is on the order of 10 −1 − 10 −2 mm a −1 (Fig. 11), while glacial erosion rates are most often on the order of 1 − 10 mm a −1 (e.g.Gurnell et al., 1996;Hallet et al., 1996;Koppes and Montgomery, 2009;Riihimaki et al., 2005).Bedrock erosion by abrasion from sediment-bearing subglacial water appears negligible compared to reported erosion rates in proglacial areas.Our results corroborate the assumption that subglacial meltwater efficiently removes sediment from the bed and we postulate that this flushing action is necessary for glacial abrasion and quarrying to access an exposed bed and remain efficient.4) through comparison of the following quantities calculated for one model year.
(a) Erosion rate averaged over the whole glacier bed ( E tot dtdx).(b) Apparent erosion rate calculated as the volume of sediment that is transported across the last grid node, i.e. terminus, ( min(q s W ch , Q tc )dt) averaged over the glacier area.This quantity corresponds to what one would measure as the sediment flux in a proglacial stream.(c) Maximum incision depth (max ėtot dt ).Simulations are ranked by averaged erosion rate and the colours represent different simulation suites: black for "reference", blue for "water input", purple for "ice geometry", red for "drainage width", and orange for "sediment supply" (Table 4).

Can ordinary seasonal melt processes lead to subglacial bedrock channel incision?
We find maximum modelled vertical bedrock incision ranging from ∼ 50 to ∼ 200 mm a −1 (Fig. 11c).Assuming that over a period of 20 years climate is relatively steady and the bedrock does not change significantly, the location of moulins would remain relatively fixed laterally and so would the channel paths (Gulley et al., 2012).Using the lowest incision rate (T_SSP/4), an N-channel almost a metre deep and a few metres wide could be carved near an ice sheet margin in 20 years.A similar N-channel would be carved in only five years assuming the largest incision rate (T_SSPOPT).Landforms created by former continental ice sheets indicate that subglacial waterways can occupy persistent paths throughout a deglaciation.Eskers deposited by the retreating Laurentide ice sheet can be traced for up to several hundred kilometres and show a dendritic pattern almost as far upstream as the former divide (e.g.Storrar et al., 2014).Some tunnel valleys show several cut-and-fill structures suggesting different excavation events; moreover, tunnel valleys carved during different glaciations tend to follow the same paths (e.g.Jørgensen and Sandersen, 2006).Eskers also commonly lie inside tunnel valleys (e.g.Jørgensen and Sandersen, 2006;Kehew et al., 2012).Assuming an incision rate of 100 mm a −1 (e.g.Fig. 7), a simple volume calculation suggests that it would take about 15 000 years to carve a 30 m deep and 100 m wide V-shaped tunnel valley, similar to the dimensions of tunnel valleys observed in Ireland (Knight, 2003).
In the context of an alpine glacier, valley geometry tends to focus subglacial water flow paths toward the thalweg.Assuming that a glacier occupies topography strongly imprinted by fluvial processes, erosion by subglacial water flow may tend to preserve if not enhance the pre-existing fluvial features along the valley centreline.For an alpine glacier eroding its bed at a pace of 2 mm a −1 (e.g.Hallet et al., 1996;Riihimaki et al., 2005;Koppes and Montgomery, 2009), in the case of simulation T_2008 (Table 3, apparent erosion of ∼ 2 mm a −1 ; ice geometry comparable to that of a large valley glacier), the maximum incision depth in one year is ∼ 125 mm (Fig. 11c).The relief of a canyon with the maximum width of the N-channel (∼ 4.5 m, for T_2008) would increase by more than ∼ 120 mm a −1 (rate of vertical bedrock incision minus rate of surrounding glacial erosion).If the canyon were 5 times as wide (∼ 22.5 m), the maximum rate of relief increase would still be ∼ 24 mm a −1 , about twice the measured incision rates in a metres-wide gorge in the French Western Alps (e.g.Valla et al., 2010), highlighting the erosional power of localized subglacial meltwater action.
5.5 What are the implications of the water flow regime in channels for hysteresis and sediment transport?
We calculate the direction of daily hysteresis between modelled transport stage (Fig. 12a) and water discharge (Fig. 12b) at four locations within the last 10 km of the glacier profile (Fig. 12c) for simulation T_2008 (Table 4; Fig 6).Overall, hysteresis is dominated by clockwise events, with anticlockwise events only occurring during the second half of the melt season.Clockwise hysteresis correlates well with the rising limb of multi-day water discharge and transport stage peaks, while anticlockwise hysteresis correlates with the falling limb, particularly at km 46 (Fig. 12).During the rising limb of a multi-day melt event, changes in channel size are dominated by enlargement; the pressure in the channel therefore peaks before the discharge, as do the averaged water flow velocity and transport stage.During the falling limb of the melt event, if the channelized drainage is relatively well established, changes in channel size are dominated by closure, and the pressure peak can occur after the peak in discharge.In the case of a proglacial stream carrying a sediment load smaller than its transport capacity, peaks in transport  stage would act as mobilizing events propagating sediment pulses downstream.We therefore surmise that the direction of hysteresis in sediment transport and discharge is not necessarily linked to changes in sediment supply conditions or the tapping of new sediment sources, but may be the result of changes in subglacial sediment mobilization in the vicinity of the glacier terminus.
In coarse-bedded streams, grain hiding has a significant effect on sediment transport (e.g.Yager et al., 2012;Schein-gross et al., 2013), as mobile grains can be trapped behind larger immobile particles.We calculate the maximum particle diameter for which movement would be initiated in simulation T_2008 (Fig. 13) and find that boulders of up to 70 cm in diameter can be transported within the last 10 km of the profile and would correspond to flood-like conditions in rivers.The analogy to river systems might have influenced the interpretation of glacial deposits such as eskers, where the presence of boulders or lack of fines is often used to infer emplacement during flood events (e.g.Brennand, 1994;Burke et al., 2012).
For simulation T_2008, we find that transport stage exhibits a sharp decrease close to the terminus (Figs.5a, 6a and 6c) which leads to a correspondingly sharp decrease in the size of particles transported (Fig. 13) and could lead to a bottleneck in sediment transport.This bottleneck effect could lead to the deposition of sediment, filling the channel toward the end of the melt season.A similar process, but operational over a longer timescale, would be consistent with time-transgressive deposition of eskers near the mouths of R-channels beneath retreating ice margins (e.g.Brennand, 1994;Burke et al., 2012).

Conclusions
This study is the first attempt to quantify bedrock erosion rates by transient subglacial water flow with a numerical model.We implement a 1-D model of subglacial drainage in which a network of cavities and R-channels interact.We compute the shear stress exerted on the bed and the resulting bedrock erosion by abrasion (saltation erosion, after Sklar and Dietrich, 2004, and total load erosion after Lamb et al., 2008).Because of the large calculated transport stage we argue that, in the case of subglacial meltwater erosion, it is probably more appropriate to use the TLEM than the SEM.Assuming that a significant amount of meltwater is produced and reaches the bed, the main drivers of subglacial water erosion that we isolate are the rate of sediment supply, particularly the relative sediment supply, and ice geometry.
From this exercise, we conclude the following: 1. Bedrock erosion and transport stage in the subglacial drainage system do not scale directly with water discharge.Instead, transport stage and erosion are related to the hydraulic potential gradient and hence a combination of water discharge, ice-surface slope, and channel (or cavity) cross-sectional area.In our simulations, this combination of discharge, slope, and channel crosssectional area leads to a drop in transport stage close to the terminus as water pressure approaches atmospheric.
2. Erosion rates due to the action of subglacial water flow averaged over the whole glacier bed are negligible compared to the rates of glacial erosion necessary to produce the sediment supply rates we impose.
F. Beaud et al.: Erosion by subglacial water flow 3.In our transient simulations, a bedrock channel a few to several decimetres in depth could be carved over a single melt season as erosion is concentrated at the base of R-channels.
4. The vertical incision rates we calculate are a few to several times larger than published rates of fluvial incision in gorges.Therefore, this mechanism may explain the gradual excavation of tunnel valleys in bedrock and the preservation or even initiation of inner gorges.
Though we have demonstrated the potential for subglacial water flow to incise bedrock on seasonal timescales, sitespecific and quantitative assessments of its importance will require more realistic 2-D hydrology models (e.g.Hewitt, 2013;Werder et al., 2013) and simulations over timescales of glacial advance and retreat.

Figure 1 .
Figure 1.(a) Different ice geometries considered (see Tables3 and 4) and (b) corresponding ice-surface slopes.The geometries REF, IS1300, and IS700 are parabolas such that the surface elevation is given by z s (x) = z s,max √x, where z s,max is the maximum ice-surface elevation and x ∈ [1, 0].In order to obtain the steeper front in STP we use a cubic root instead of the square root: z s (x) = z s,max x 1/3 .The wedge-like geometry WDG is defined by a straight line with the same elevation change as REF.

Figure 2 .Figure 3 .
Figure 2. Relationship between the hydraulic and erosion variables in a steady-state R-channel with constant discharge (Table3, S_MOULIN; Sect.4.1.1).Components of TLEM and SEM are also compared.(a) R-channel discharge, Q ch , and cross-sectional area S; (b) velocity, u, and gradient in hydraulic potential, ∇φ ch ; (c) sediment transport capacity, q tc , prescribed sediment supply rate, q s , and fraction of bed exposed, F e ; (d) near-bed sediment concentration, c b (TLEM), and impact rate, I r (SEM); (e) impact velocity in the TLEM, w i,eff , and the SEM, w si ; and (f) erosion rate calculated with the TLEM, ėtot , and the SEM, ėsalt .

Figure 4 .
Figure 4. Influence of sediment supply rate q s (S_SSP, Table3) on total erosion with the TLEM and SEM for the same hydrology as the reference simulation (S_REF, Table3).The legend applies to all panels.(a) Relative sediment supply (q s /q tc ), (b) total erosion rate computed with the TLEM (E tot ), and (c) total erosion rate computed with the SEM (E salt ).Note the logarithmic scale on the y axis of (b) and (c).

Figure 8 .
Figure 8.Comparison of time-integrated erosion per unit width ( ėtot dt) and total erosion ( E tot dt) for different ice geometries.The first 20 km of the profile are not shown because the bed is alluviated and erosion rates are negligible.

Figure 9 .
Figure 9. Influence of drainage catchment width (W ) on timeintegrated erosion.(a) Time-integrated erosion per unit width ( ėtot dt) and total erosion ( E tot dt) for an individual R-channel.(b) Time-integrated total erosion ( E tot dt) summed over all Rchannels in each simulation.The first 20 km of the profile are not shown because the bed is alluviated and erosion rates are negligible.
Figure10.Time-integrated erosion per unit width ( ėtot dt) and time-integrated total erosion ( E tot dt) as a function of sediment supply rate (q s ) (Table4).The first 20 km of the profile are not shown because the bed is alluviated and erosion rates are negligible.

FFigure 11 .
Figure 11.Synthesis of transient simulations (Table4) through comparison of the following quantities calculated for one model year.(a) Erosion rate averaged over the whole glacier bed ( E tot dtdx).(b) Apparent erosion rate calculated as the volume of sediment that is transported across the last grid node, i.e. terminus, ( min(q s W ch , Q tc )dt) averaged over the glacier area.This quantity corresponds to what one would measure as the sediment flux in a proglacial stream.(c) Maximum incision depth (max ėtot dt ).Simulations are ranked by averaged erosion rate and the colours represent different simulation suites: black for "reference", blue for "water input", purple for "ice geometry", red for "drainage width", and orange for "sediment supply" (Table4).

Figure 12 .
Figure 12.Hysteresis between water discharge and transport stage for simulation T_2008 (see Table4and Fig.6).Time series at four locations(km 40, 43, 46, and 50)  distributed over the last 10 km of the glacier profile of (a) transport stage (τ * /τ * c ) and normalized water input ( ḃca / ḃss max ) at the terminus (light grey); (b) discharge in the channel, Q ch ; and (c) calculated direction of the daily hysteresis when transport stage is plotted against water discharge.The hysteresis is clockwise when transport stage peaks before discharge over a daily cycle.Undefined events represent days where the fluctuations in transport stage or discharge are either simultaneous or not strong enough to produce hysteresis.

Table 1 .
Summary of hydrological model parameters (see Supplement for extended list).
Comparison of time-integrated erosion per unit width ( ėtot dt) and total erosion ( E tot dt) for different water inputs.The first 20 km of the profile are not shown because the bed is alluviated and erosion rates are negligible.remains highest at 42.5 km because creep closure prevents R-channels from becoming too large and water discharge is high enough to maintain high velocities.At km 20 and 27.5 it takes about 60 days for the first peak in transport stage to occur.The subsequent peaks at these two locations (around day 90, 105, and 129) lag high melt periods even further.