Modeling debris-covered glaciers: response to steady debris deposition

. Debris-covered glaciers are common in rapidly eroding alpine landscapes. When thicker than a few centimeters, surface debris suppresses melt rates. If continuous debris cover is present, ablation rates can be signiﬁcantly reduced leading to increases in glacier length. In order to quantify feedbacks in the debris–glacier–climate system, we developed a 2-D long-valley numerical glacier model that includes englacial and supraglacial debris advection. We ran 120 simulations on a linear bed proﬁle in which a hypothetical steady state debris-free glacier responds to a step increase of surface debris deposition. Simulated glaciers advance to steady states in which ice accumulation equals ice ablation, and debris input equals debris loss from the glacier terminus. Our model and parameter selections can produce 2-fold increases in glacier length. Debris ﬂux onto the glacier and the relationship between debris thickness and melt rate strongly control glacier length. Debris deposited near the equilibrium-line altitude, where ice discharge is high, results in the greatest glacier extension when other debris-related variables are held constant. Debris deposited near the equilibrium-line altitude re-emerges high in the ablation zone and therefore impacts melt rate over a greater fraction of the glacier surface. Continuous debris cover reduces ice discharge gradients, ice thickness gradients, and velocity gradients relative to initial debris-free glaciers. Debris-forced glacier extension decreases the ratio of accumulation zone to total glacier area (AAR). Our simulations reproduce the “general trends” between debris cover, AARs, and glacier surface velocity patterns from modern debris-covered glaciers. We provide a quantitative, theoretical foundation to interpret the effect of debris cover on the moraine record, and to assess the effects of climate change on debris-covered glaciers.


Introduction
Glaciers erode landscapes directly by subglacial quarrying and abrasion, and indirectly by steepening hillslopes above glaciers. Oversteepened hillslopes can deliver loose rock (debris) onto glacier surfaces (Benn and Evans, 2010). Steep hillslopes and high hillslope erosion rates in alpine settings therefore tend to correspond with the occurrence of debriscovered glaciers (e.g., the Himalaya and the Alaska Range). We refer to a debris-covered glacier as any glacier with continuous debris cover across the full glacier width over a portion of the glacier (after Kirkbride, 2011).
Debris cover more than a few centimeters thick damps the ablation of underlying ice (e.g., Østrem, 1959;Shroder et al., 2000;Owen et al., 2003). If debris supply to a glacier surface is high, mass balance profiles can be greatly altered, leading to increases in glacier volume and length (e.g., Konrad and Humphrey, 2000;Scherler et al., 2011a, b;Fig. 1). Thick debris cover on glaciers can also lead to low accumulation-area ratios (AARs; Scherler et al., 2011b). Estimates of past climate change will therefore be exaggerated if typical AARs are assumed when reconstructing past climate from moraines deposited by debris-covered glaciers.
Debris-covered glacier termini exhibit a wide range of responses to climate change (Scherler et al., 2011a). While almost all Himalayan debris-free glaciers are retreating, Himalayan debris-covered glacier termini are not respond- Figure 1. (a) Schematic of the debris-glacier system. Debris deposited on or emerging in the ablation zone reduces ablation rates (above the critical debris thickness) leading to the reduction in gradients of ice discharge and the lengthening of glaciers. (b) Schematic of the coupled debris-glacier model. Debris deposited on the glacier is either advected through the glacier and/or advected down the glacier surface. Englacial debris is advected using 2-D rectangular grid and coordinate transform. Ice physics and supraglacial debris advection is treated on a 1-D grid.
ing coherently to climate change despite a strong trend toward negative mass balance and surface lowering (e.g., Bolch et al., 2011Bolch et al., , 2012Benn et al., 2012;Kääb et al., 2012). Some Himalayan debris-covered glacier termini are advancing, others are stationary, and yet others are retreating (e.g., Raper and Braithwaite, 2006;Scherler et al., 2011a;Bolch et al., 2012). This discrepancy between debris-covered glacier mass balance and terminal response highlights the pressing need to understand the sometimes counterintuitive effects of debris on glacier response.
The direct effect of debris on glaciers is difficult to isolate on modern glaciers. In situ documentation of debris-covered glacier mass loss is made difficult by non-uniform debris thicknesses and the presence of scattered ice cliffs, surface ponds, and proglacial lakes. As a result complete summer balances from debris-covered glaciers are sparse (WGMS, 2008). Measurements of englacial debris concentrations and distribution are yet more difficult to obtain (e.g., Kirkbride and Deline, 2013). In addition, exploration of century-scale response of debris-covered glaciers is limited by short satellite and observational periods (Bolch et al., 2011). Logistical realities therefore limit our ability to constrain feedbacks between debris deposition rates, the englacial environment, the supraglacial environment, and ice dynamics.
While logistics limit our ability to directly observe some feedbacks, many of the most provocative conclusions relating debris and glacier response are based on remotely sensed data. Scherler et al. (2011b) provided an extensive inventory of remotely sensed velocity and debris coverage data from 287 glaciers in High Asia. They inferred several general patterns from these debris-covered glaciers: (1) hillslope debris flux onto glaciers correlates with the percentage of debris cover on glaciers; (2) debris-covered glacier AARs tend to be smaller than debris-free glaciers; and (3) surface debris perturbs velocity distributions on valley glaciers by shifting maximum glacier velocities up glacier, away from the terminus. These inferences highlight the effect of thick debris cover on valley glaciers and serve as targets for models of debris-covered glaciers.
Numerical models can help quantify feedbacks within the debris-glacier-climate system (e.g., Konrad and Humphrey, 2000). Debris-covered glacier models have been used to explore the response of valley glaciers to (1) the steady input of debris (Konrad and Humphrey, 2000); (2) one-time landslide deposition of debris on glaciers (Vacco et al., 2010;Menounos et al., 2013); and (3) climate change (Naito et al., 2000;Banerjee and Shankar, 2013;Rowan et al., 2015). Konrad and Humphrey (2000) used a twodimensional (2-D; long-valley-vertical) model with a constant surface slope to explore debris-covered glacier dynamics. In their model, debris was deposited on the glacier surface below the equilibrium-line altitude (ELA) and was then advected along the glacier surface. With high debris fluxes, simulated glaciers formed several-meter-thick debris covers, which reduced sub-debris melt toward zero, and resulted in glaciers that never reached steady state. Numerical models have also shown that large landslides onto glaciers can lead to multiple-kilometer advances of the terminus (Vacco et al., 2010;Menuounos et al., 2013). Debris-covered glacier retreat response timescales have also been explored with a simplified debris-covered glacier model (Banerjee and Shankar, 2013). Rowan et al. (2015) used a numerical model to forecast the response of the debris-covered Khumbu glacier, Nepal to climate change. But owing to the complexity of the debris-glacier-climate system, it can be difficult to diagnose the effects of different processes on observable glacier responses. For example, both increased debris delivery to a glacier and a cooling climate could lead to glacier advances (e.g., Vacco et al., 2010;Menuounos et al., 2013). What approaches could we use to address these sorts of conundrums within the debris-glacier-climate system?
Here we attempt to improve our understanding of the debris-glacier-climate system (and subsequently better project future glacier change) by isolating how debris effects glacier response, while holding climate steady. While significant effort has focused on glacier-climate interaction, less research has focused on isolating the effect of debris on glacier length (e.g., Konrad and Humphrey, 2000), and other basic measures of glacier response (e.g., change in glacier surface velocity due to debris deposition on the glacier). We explore debris-glacier interactions by isolating the role of debris in governing basic glacier dynamics and glacier length.
We use a simple glacier model to simulate hypothetical debris-covered glaciers. This new framework allows us to isolate the effects of debris on glacier response by controlling the potentially conflating effects of a variable bed, variable glacier width, or a temporally variable climate. To isolate the effect of debris, we start each simulation with a steady state debris-free (ssdf) glacier and impose a step change increase in debris deposition rate while holding climate steady. In many debris-covered glacier systems, debris is deposited in the accumulation zone, advected through the glacier, and emerges in the ablation zone (e.g., Boulton and Eyles, 1979;Owen and Derbyshire, 1989;Benn and Owen, 2002;Benn et al., 2012). Our new transient 2-D numerical model (x, z) links debris deposition, englacial debris advection, debris emergence, surface debris advection, debris-melt coupling, debris removal from the glacier terminus, and shallowice-approximation dynamics (Figs. 1 and 2). We provide a new terminus parameterization which allows for the use of steady state glacier length as a metric for comparison between simulated debris-covered glaciers. While real debriscovered glaciers may not reach steady state the concept is necessary for determining the sensitivity of debris-covered glaciers to changes in debris-related parameters. Our intent is to determine which parameters and parameterizations are most important for capturing the response of glaciers to debris input. Here, we explore the sensitivity of hypothetical debris-covered glaciers to changes in debris-input related variables (e.g., debris flux, debris deposition location, and debris deposition zone width). We also explore the sensitivity of debris-covered glaciers to different debris thickness-melt formulations. We compare our theory-based results to the "general trends" documented by Scherler et al. (2011b). By isolating the effect of debris on glaciers, this study lays the theoretical foundation for efforts exploring the more complex response of debris-covered glaciers to climate change.

Theory and numerical methods
We employ a fully transient 2-D finite difference numerical model (in down-valley and vertical, x and z) that can simulate the evolution of temperate valley glacier response to climate and debris. Forced by a time series of equilibrium-line altitudes (ELAs) and a prescribed mass balance gradient, the model calculates ice surface elevations above a longitudinal profile by solving equations for ice flux and mass conservation. The modeled longitudinal path represents the glacier centerline. A number of authors have used the shallow-iceapproximation (SIA) and basal sliding parameterizations in numerical glacier models (e.g., Nye, 1965;Budd and Jensen, 1975;Oerlemans, 1986;MacGregor et al., 2000;Leysinger and Gudmundsson, 2004;Kessler et al., 2006). We employ a similar approach, but add a longitudinal stress coupling parameterization (Marshall et al., 2005). The model is efficient, allowing wide exploration of parameter space in simulations over thousands of years.

Conservation of ice mass
Mass conservation is at the core of the ice physics model. Assuming uniform ice density, and ignoring variations in the width of the glacier, ice conservation requires that where x is the distance along the glacier flowline, H is the local ice thickness,ḃ is the local specific balance, and Q [in m 3 m −1 yr −1 ] is the specific volume discharge of ice. This requires a prescribed mass balance field, and a prescription of the ice physics governing ice discharge. We use a simple mass balance scheme that limits the number of parameters while honoring the essence of glacier surface mass balance. We combine surface accumulation and ablation into a single thresholded net mass balance profile as a function of elevation, z: where dḃ z dz is the mass balance gradient with elevation, Z ice is the ice surface elevation andḃ max z is a maximum mass balance that accounts for the depletion of moisture available for precipitation at higher elevations. The annual surface mass balance of ice in the absence of debris is held steady for all simulations to isolate the effects of debris from those of climate change on glacier response.

Annual surface mass balance: effect of supraglacial debris
Sub-debris melt rate decreases rapidly with increasing debris thickness (e.g., Østrem, 1959;Nicholson and Benn, 2006). For debris layers thinner than a critical thickness (∼ 2 cm), surface debris can increase melt rates relative to bare ice. For debris thicknesses greater than ∼ 2 cm, debris suppresses sub-debris melt rates relative to bare ice (e.g., Nicholson and Benn, 2006;Fig. 3). We assume that heat is transferred through the debris layer by conduction. Sub-debris melt should therefore vary inversely with debris thickness (i.e., be hyperbolic) as conduction is governed by the temperature gradient ∼ (T s − T ice )/h debris (e.g., Nicholson and Benn, 2006). Here, T ice = 0. We neglect the melt-amplifying effects of very thin debris for simplicity and represent the damping of sub-debris melt rates with where h * is a characteristic length scale and k and φ are thermal conductivity and porosity of debris cover, ρ i and L the density and latent heat of fusion of ice, T s the average debris surface temperature, T a the average screen-level air temperature, and f pdd is a positive degree day factor relating air temperature and the bare ice melt rate (e.g., Mihalcea et al., 2006). In this formulation, subdebris melt rates approach bare-ice melt rates as debris thins (h debris h * ), and asymptotes toward a hyperbolic dependence on debris thickness as debris thickens (h debris h * ).
For comparison, we also show the most likely exponential fit to the data (Fig. 3). The exponential curve fit declines more rapidly than the hyperbolic fit (e.g., Konrad and Humphrey, 2000;Hagg et al., 2008). We neglect the effects of surface streams, thermokarst, and ice cliffs that can lead to complex local topography and melt rates within debris covers (e.g., Reid and Brock, 2014;Anderson, 2014).

Ice dynamics
Ice is transferred down valley by internal ice deformation and by basal motion. The ice discharge down glacier is in which H is the local ice thickness and u is the depthaveraged bed parallel velocity that results from the sum of the ice deformation velocity and basal motion. The SIA reduces the momentum balance equations to expressions for vertical shear stress as a function of the local ice surface slope and ice thickness. The depth-averaged horizontal velocity due to internal deformation is where ρ i the density of ice, g the acceleration due to gravity, α the local ice surface slope, H the local ice thickness, τ bx is the local basal shear stress, A is the creep parameter, and n is the flow law exponent (assumed to be 3). We assume that all ice is temperate, and A is therefore taken to be 24×10 −25 [in Pa −3 s −1 ] (Cuffey and Paterson, 2010). In addition to internal deformation, temperate glaciers transfer mass via basal slip due to ice sliding over the bed and deformation of the bed itself. We assume that all basal slip is accomplished by sliding over bedrock, and follow the formulation of Kessler et al. (2006): in which u c is a typical sliding velocity, and τ c is the gravitational driving stress that gives rise to the typical sliding velocity. This sliding parameterization is not as sensitive to high τ b values as many other sliding laws (e.g., Cuffey and Paterson, 2010), and provides a more conservative estimate of sliding velocities when τ b > τ c (Kessler et al., 2006). We have also modified the SIA equations by including a parameterization of longitudinal stress coupling (after Marshall et al., 2005) and a shapefactor, f , that represents the effect of valley wall drag. The longitudinal coupling scheme modifies τ bx to where the effective viscosity, η = 1 2 [Aτ n−1 E ] −1 and u is the vertically averaged ice velocity. In the shallow ice approximation, τ E , the effective stress, is approximated by the local τ bx (after Cuffey and Paterson, 2010). We take f = 0.75 to approximate the effects of sidewall drag from a parabolic valley cross-section with a half-width 3 times the ice thickness (Cuffey and Paterson, 2010).

Ice velocity structure within the glacier
Horizontal and vertical velocity fields must be resolved within the glacier in order to advect englacial debris. We start by defining the horizontal velocity field within the glacier, and then employ continuity in an incompressible medium to calculate the associated vertical velocities. The u(z) profile shape may be obtained from the analytic solution to flow of ice in a uniform channel with Glen's flow law rheology: where ζ is the non-dimensional height z/H above the bed, and F = u(z) u def is the ratio of horizontal speed to mean deformation speed. The full horizontal velocity field is then characterized by where u coupling is the vertically integrated velocity effect due to longitudinal stress coupling and is determined by subtracting the original Eq. (6) from Eq. (6) modified by Eq. (8).
Vertical and horizontal velocity fields (w(x, z) and u(x, z)) are related through the continuity equation for an incompressible fluid, which in two dimensions (x, z) is We then solve for the vertical velocity in each cell within each column by integrating vertically, as follows: employing the boundary condition that w = 0 at z = 0 (i.e., we assume no basal melt). In steady state, vertical velocities, w, at the glacier surface must be equal in magnitude and opposite in sign to the surface mass balance field, and are therefore directed downward at the ice surface in the accumulation zone, and upward in the ablation zone.

Debris deposition
Debris can be entrained in the glacier at the upper glacier surface, at the glacier bed, or along glacier sidewalls (e.g., Hambrey et al., 2008). Supraglacial debris deposition largely occurs by mass wasting from hillslopes above glaciers, while sub-glacial/sidewall debris entrainment occurs through regelation and net freeze-on. Basal debris emergence at the glacier surface is typically limited to the glacier toe and likely plays a minor role in the formation of extensive debris covers (Benn and Evans, 2010). We focus on debris sourced from valley head and side walls. Headwall erosion rates are better constrained than subglacial entrainment rates and mass wasting from head and sidewalls is the primary process of debris delivery onto many valley glaciers (Messerli and Zurbuchen, 1968;Humlum, 2000 (Benn and Evans, 2010). For simplicity, we neglect englacial thrusting and ice-stream interaction moraines (medial moraines associated with tributary junctions; see Eyles and Rogerson, 1978;Anderson, 2000;Benn and Evans, 2010). These cases can be treated in subsequent modeling that incorporates the 2-D planform complexities of valley glaciers. Debris delivery to glacier surfaces can vary considerably from glacier to glacier, depending on glacier topology and above-glacier topography (e.g., Deline, 2009). We capture this complexity using four variables: the total debris flux to the glacier surface (ḋ flux , [in m 3 m −1 yr −1 )], the debris deposition rate (ḋ [in mm yr −1 )], the debris deposition zone width (d width [in m]), and the debris deposition location (d loc ). In the model,ḋ flux is representative of the integrated effects oḟ d and d width .
Rock type, slope, and fracture density are significant factors determining hillslope erosion rates and therefore also control the debris deposition rate,ḋ (e.g., Stock and Montgomery, 1999;Molnar et al., 2007). In the model,ḋ, is allowed to vary from 1 to 8 mm yr −1 and is steady within each simulation (Fig. 1b). Debris deposition rate depends on a number of site-specific variables: where f funneling is a dimensionless factor capturing the effect of topographic funneling on debris deposition, f hillslope is the percentage of the headwall that is exposed bedrock, e is the hillslope backwearing rate in m yr −1 , H wall is the height of the headwall, and θ is the headwall slope. The deposition rates explored in this study are appropriate for typical headwall erosion rates (typically ranging between 0.5 and 2 mm yr −1 ), headwall heights, and headwall slopes for highrelief mountain environments (e.g., Heimsath and McGlynn, 2008;Ouimet et al., 2009;Scherler et al., 2011;Ward and Anderson, 2011). d width defines the down-valley width of the deposition zone, the zone over which the debris is spread on the glacier surface (we employ a base width of 400 m; Table 1; Fig. 1b). Debris is deposited on glaciers at locations where hillslope erosion processes are connected to the glacier surface. This requires high-relief topography above the glacier to provide the energy necessary to move the debris onto the glacier. In the model, we control the down-valley debris deposition location with the variable d loc , which we allow to vary from near the headwall to near the glacier terminus. d loc defines the up glacier end of the debris deposition zone.

Incorporation and advection of englacial debris
Debris deposited in the ablation zone is advected along the glacier surface, whereas debris deposited in the accumulation zone moves downward with the ice and is therefore incorporated into the glacier. The near-surface debris concentration in the accumulation zone is defined as where m z is the number of vertical slices the englacial advection scheme is divided into (H /m z being the thickness of the slices) and dt is the model time interval. C 0 is therefore the mass concentration of debris in the surface-bounding cell.
Once embedded in the glacier, debris is advected through the glacier following englacial flowpaths ( ∂C ∂t = − ∂(uC) ∂x − ∂(wC) ∂z ). Taking an Eulerian point of view, the time rate of change of concentration of debris within an ice cell (in our model) is where h ζ is the cell height in a given ice column (h ζ = H m z ). The first and second terms represent changes in C due to advection in the vertical and the horizontal directions, respectively. The third term on the right hand side represents the rate of change of C due to vertical ice strain from the thinning or thickening of the glacier through time. Note that if the strain rate is negative, signifying vertical thinning of an ice column, debris concentration in a cell will increase. The fourth term represents the rate of change of C due to the longitudinal changes in glacier thickness. This term accounts for the fact that cells from one column to the next are not the same volume.

Advection of debris on the glacier surface and steady states
We track both the melt-out of englacial debris and the advection of supraglacial debris on the glacier surface. The rate of change of debris thickness on the glacier surface is captured by where h debris is the debris thickness, ρ rock is the density of the rock, φ is the porosity of supraglacial debris, and u surf is the surface velocity of the glacier (after Konrad and Humphrey, 2000;Naito et al., 2000;Vacco et al., 2010). The first term on the right represents the addition of debris to the surface from melt of debris-laden ice. The second term represents the advection of debris down glacier. Debris is transported off glacier by the wasting of debris down the terminal slope or by the backwasting of terminal ice cliffs (Konrad and Humphrey, 2000;Appendix A and B). In the model we implement a triangular terminus wedge parameterization (after Budd and Jenssen, 1975; see Appendix A). The change of surface debris thickness with time on the terminal wedge is whereḋ term flux is the debris flux into the foreland from the terminus wedge [in m 3 m −1 yr −1 ] and dx term is the surface length of the terminal wedge. We useḋ term flux =ḃ term z h term debris . Varying this parameterization has a minor effect on glacier length, but can have a considerable effect on the temporal evolution of the glacier asḋ flux must equalḋ term flux for a simulated glacier to reach steady state (Appendix A). We explore the choice and effect of this parameterization in Appendix B.

Implementation and numerics
We now outline the order of calculations in the model. First,ḃ z and b are calculated based upon elevation and debris thickness. Next, we use a second-order Runge-Kutta centered difference scheme to evolve H (x, t), followed by the implementation of an iterative "upstream" debris advection scheme following Smolarkiewicz, 1983. This iterative scheme imposes a two-step anti-diffusion correction algorithm to the advection scheme, which greatly reduces numerical diffusion (Smolarkiewicz, 1983). We test advection scheme stability using the Courant-Friedrichs-Lewy (CFL) condition, which ensures that mass is not advected beyond adjacent cells in a single timestep. We implement a terminus wedge parameterization that allows simulated glaciers to advance to steady state (Appendix A and B). The time step, dt, for ice-physics and debris advection is 0.01 years. All ice columns are segmented into m z heights (i.e., ζ = 0 : (1/m z ) : 1); in all results below we use m z = 20 (Fig. 1b). We impose a no-flux boundary at the upper end of the glacier. While our simulations are hypothetical we select the base model parameters to loosely represent the ablation zones of debris-covered glaciers in the Khumbu region of Nepal. There is a wealth of debris-covered glacier research from this region, which assures that our parameter choices in the range of observed values (e.g., Kayastha et al., 2000;Bolch et al., 2011;Benn et al., 2012;Shea et al., 2015). Base simulations are run on a linear glacier bed with a basal slope of 8 % and a maximum bed elevation of 5200 m (Scherler, 2014). This simple bed geometry is used to ensure that our results to do not conflate the effects of bed topography with the effects of debris. We use a dḃ dz = 0.0075 yr −1 , which is capped at 2 m yr −1 based on data from debris-free glaciers in the Khumbu region (Mera and Pokalde glaciers: after Wagnon et al., 2013). Our parameter exploration below shows that our conclusions are not influenced by our choice of base parameters from the ablation zones of debris-covered glaciers in the Khumbu region. All simulations start with an 8.7 km long steady state debris-free (ssdf) glacier with a steady ELA at 5000 m (L ssdf = 8.7 km). In each simulation a step change increase in debris deposition rate is imposed at t = 100 years. The base parameter set usesḋ flux = 3.2 m 3 m −1 yr −1 ,ḋ = 8 mm yr −1 ,ḋ width of 400 m, and the location of debris input, d loc , is 42 % of the distance between the headwall and the length of the debris-free glacier, L ssdf .
www.the-cryosphere.net/10/1105/2016/ The Cryosphere, 10, 1105-1124, 2016 We first demonstrate the transfer of debris between model components and demonstrate debris-covered steady state. We then explore the differences between the steady state debrisfree (ssdf) glacier and debris-covered glaciers and explore relative importance ofḋ, d width , d loc ,ḋ flux , andḋ term flux on glacier length. The effect ofḋ term flux on the length and time evolution of the model is explored in Appendix B (see Fig. B1). We then test the sensitivity of the model to changes in h * and φ. Last, we compare our hypothetical simulations to "general trends" observed from real debris-covered glaciers.

Demonstration of debris-covered glacier steady state and conservation of debris
In order to compare steady state glacier lengths between simulations with differentḋ flux we track debris through the model. At any time in the simulation, the total debris mass that has been deposited on the simulated glacier must equal the total debris mass in the model: where M input is the total rock mass deposited on the glacier and accumulated over time, M englacial is the total englacial debris mass, M surface is the total debris mass on the glacier surface, and M foreland is the total mass deposited in the proglacial environment. We use the base parameter set simulation to highlight the transfer of debris mass through the system (Fig. 4). Because debris is deposited in the accumulation zone near the ELA, in the base simulation, M englacial rapidly reaches steady state (Fig. 4). As the glacier extends, M surface continues to increase at a declining rate as more surface debris is transferred into the foreland. The glacier reaches steady state when the glacier length, M surface , and M englacial are steady and the rate of change of M foreland is equal the rate of debris input to the glacier. Each model simulation presented conserves greater than 99 % of debris mass.

Comparison of modeled debris-free and debris-covered glaciers
We first highlight differences in length, and the patterns of ice discharge, Q, ice thickness, H , and surface speed, u surf , between the ssdf glacier and its steady state debris-covered counterpart, using the base parameter set (Fig. 4). In this baseline case the steady state debris-perturbed glacier length is 175 % of L ssdf (Fig. 5).
The debris thickness, h debris , increases down glacier from the site of debris emergence, x int , except near the terminal wedge where theḋ term flux parameterization reduces h debris (Figs. 5-6). Down glacier from the site of debris emergence, x int , gradients of Q, H , and u surf are reduced relative to the debris-free glacier ( Fig. 6b and d). Debris-free patterns of Figure 4. Debris mass vs. time. The englacial debris mass reaches steady state rapidly because debris is deposited near the ELA and englacial advection paths are short. As debris emerges in the ablation zone M surface increases nearly at the rate of debris input to the glacier. As the glacier nears a steady length the debris mass transferred to the glacier foreland increases. The glacier reaches steady state whenḋ flux =ḋ term flux and the glacier length is steady (see Appendix A).  Q and u surf are convex up near the glacier terminus, while Q and u surf from debris-covered termini are concave upward. The lowest gradients in Q, H , and u surf occur near the glacier terminus where h debris is thickest (excluding the terminal slope; Fig. 6).

Effect of debris input location
Debris input location (d loc ) controls the englacial debris path. Debris deposited near the headwall is advected more deeply into the glacier than debris deposited near the ELA. Debris deposited near the ELA follows a shallow, short englacial path (Fig. 5). The original width of the debris band deposited in the accumulation zone, is reduced down glacier and then widens again near the surface in the ablation zone (Fig. 5). The debris band initially narrows due to the longitudinal straining of ice (Hooke and Hudleston, 1978;Cuffey and Paterson, 2010;Fig. 5a) and then widens due to feedbacks between the surface debris and ice dynamics.
In order to show the effects of d loc on basic glacier properties (glacier length, Q, H , and u surf ), we highlight three simulations where we vary d loc and hold all other debrisrelated parameters constant (ḋ flux = 3.2 m 3 m −1 yr −1 ,ḋ = 8 mm yr −1 , andḋ width = 400 m) . d loc is varied from near the top of the glacier to near the debris-free glacier toe (Figs. 5 and 6).
When debris is deposited or emerges where Q is large (near the ELA), glacier extension is greater than when debris is deposited/emerges where Q is small (near the headwall or the debris-free glacier terminus). If debris is deposited or emerges where Q free /Q max nears 1 glacier extension will be largest for a given glacier (Q free refers to ice discharge from the ssdf glacier and Q max is the maximum Q free before debris is added to the glacier). Where Q free /Q max nears 0 glacier extension will be small.
We ran an additional 33 simulations (36 total) in which we varyḋ flux and d loc (Fig. 7). Varying the debris deposition location while holding the debris flux constant results in a maximum of a 40 % difference (for these 36 simulations; Table 2) in the resulting steady-state debris-covered glacier length. The importance of d loc on glacier length increases with largerḋ flux (Fig. 7). The general pattern seen in Fig. 7 is insensitive to changes in other parameters. Increasingḋ flux leads to increases in the percentage of the glacier covered with debris (Fig. 8).

Effect of debris deposition rate, debris deposit width, and debris flux
Increasing either the debris deposition rate (ḋ) or the debris deposit width (d width ) leads to increases inḋ flux , but the relative importance ofḋ or d width in governing glacier response is unclear. Does debris delivered to a small portion of a glacier at a high rate lead to a different length response than debris delivered to a glacier in a wide section but at a low rate? In order to parse the effects ofḋ and d width on glacier length, we ran simulations in which we varyḋ and d width (Fig. 9b for d loc = 42%). The effect is small, varying the contribution results in a maximum of a 4 % difference in steady-state debriscovered glacier length for any given debris flux ( Fig. 9b; Table 2). In contrast, varying the debris flux,ḋ flux , results in a maximum of 80 % change in glacier length ( Fig. 9c; Table 2).

Effect of characteristic debris thickness and surface debris porosity
We explore the sensitivity of the model to changes in the characteristic debris thickness (h * ) and surface debris porosity (φ). We vary h * and φ, impose a step change increase in debris input to the ssdf glacier and compare the resulting steady state glacier lengths (Fig. 10). Simulated glacier length is highly sensitive to h * (Fig. 10). For the same debris delivery variables, the more rapidly the melt rate is damped by debris (lower h * ), the longer the steady state glacier. Steady state debris-covered glacier length varies by 110 % relative to L ssdf when h * is varied from the extremes of 0.0035 to 0.165 m (Table 2; 55 % for the 1 σ range (0.037-0.095 m)). Glacier length is not as sensitive to the choice of debris porosity, φ (Fig. 10). Varying φ between the extremes of 0.18 and 0.43 (e.g., Bozhinskiy et al., 1986; Conway and Rasmussen, 2000) leads to lengths that vary 25 % relative to L ssdf (Table 2).

Comparison with trends observed from debris-covered glaciers
Our model results show that steady, high debris fluxes onto glaciers lead to increased glacier lengths and high percentages of debris cover (Figs. 8 and 9). Remote-sensing derived measurements provide general insight into valley glacier response to debris. We compare our hypothetical results to the broad trends Scherler et al. (2011b) inferred from their inventory of 287 debris-covered glacier surface velocity patterns, AARs, and debris cover percentages. While the Scherler data set was collected from glaciers responding to persistent negative mass balance, the authors note that their inferences stand even "when excluding stagnating glaciers". This suggests that their observations represent "general trends" relating debris to glacier response (e.g., increasing debris flux leads to reduced AARs and an up glacier shift of maximum glacier surface velocities).  Scherler et al. (2011b) documented that higher debris cover percentage on glaciers correlates with steep aboveglacier hillslopes. Because hillslope erosion rates and the percentage of exposed bedrock in the headwall increase with steeper slopes, it follows that increased debris input onto a glacier should also increase both the glacier length and the percentage of the glacier covered with debris. Our hypothetical model results confirm this inference and show that -independent of parameter selection (e.g., d loc , h * , bed slope) -higher debris flux leads to higher debris cover percentages on glaciers (Figs. 8 and 11). Scherler et al. (2011b) showed that large debris cover percentages correspond with small AARs outside the typical range of 0.5-0.7 from debris-free glaciers (e.g., Meier and Post, 1979). Our modeled steady state debris free glacier has an AAR of 0.5 (due to the piecewise-linear mass balance profile and constant width valley). In our model simulations, increases in debris flux lead to increases in both steady state glacier length, and debris cover percentage independent of parameter selection (Fig. 11a). With a fixed ELA, the AAR must therefore decrease with an increased debris flux (Fig. 11a). Varying h * (using the base parameter set; Fig. 10) has a similar effect to varying debris flux ( Fig. 11c and d). Changes in the location of debris input lead to small changes in AAR but considerable changes in debris cover percentage (Fig. 11a). Scherler et al. (2011b) also showed that larger debris cover percentage correlated with lower ratios of average surface speed(u surf ) from the lower half of glaciers to the average u surf from the upper half of glaciers. Increasing the debris flux in our model leads to lower u surf in the lower half of glaciers relative to u surf in the upper half of glaciers independent of parameter selection (Fig. 11b). Changing the location of debris input, d loc , leads to small changes in the ratio of average u surf but leads to large changes in the percentage of the glacier covered with debris. This highlights that debris flux Figure 11. Comparison of our hypothetical steady state debris-cover model output with data from 287 glaciers showing broad patterns between debris and basic glacier properties (Scherler et al., 2011b). (a) The AAR compared to debris cover percentage, debris flux (ḋ flux ), and debris deposition location (d loc ). (b) The ratio of the average surface speed of the lower 50 % of the glacier and the average surface speed of the upper 50 % of the glacier vs. debris cover percentage,ḋ flux , and d loc . (c, d) Same data as (a, b), but exploring the effect of changing the bed slope and h * . The quadrangles show the area occupied by simulation results using the same parameters from (a, b) but with lower and higher bed slopes. h * results are from the parameter test where h * is varied (Fig. 10). and debris deposition location are important parameters for the specific response of a glacier to debris input.
In order to show the generality of inferences made by Scherler et al. (2011b), we also change the bed slope in our hypothetical model. Changing the linear bed slope leads to similar relationships between debris cover percentage, AAR, and surface velocity. Notable differences occur primarily when the bed slope is reduced ( Fig. 11c and d). With a reduced bed slope the initial debris-free steady state glacier is 3 times longer than the steady state debris free glacier. Even with the same hillslope debris fluxes as the simulations in Fig. 11a and b, the reduced bed slope leads to reduced asymmetry in the steady state debris-covered glacier surface velocities (Fig. 11d). With a linear mass balance profile and linear bed slope, changing the bed slope will have a similar effect to changing the mass balance gradient. The specific relationship of glacier response to debris is therefore also dependent on glacier size, bed slope, and the environmental mass balance gradient. Ultimately, our exploration shows that, independent of parameter selection (e.g., not dependent on bed slope or mass balance profile selection), our model reproduces basic patterns inferred from real debris-covered glaciers, which lends support to our model framework, while also providing quantitative, theoretical support to previous data-based observations.

Discussion
We explored the sensitivity of a new debris-covered glacier model to changes in various parameters and debris input related variables. We used a rigorous steady state glacier length definition to allow for the intercomparison of each simulation. Simulated glacier lengths are most sensitive to hillslope debris flux and the selection of the debris thickness that characterizes the decline in melt rate beneath debris (Table 2). The location of debris deposition is important but plays a secondary role in setting glacier length. The time evolution of debris-covered glacier length is highly dependent onḋ term flux , although steady state glacier length is not (Appendix B; Fig. B1). Thick debris cover on glaciers from consistent debris input, independent of climate change, tends to (1) reverse and reduce mass balance gradients; (2) extend glaciers; (3) reduce AARs; and (4) reduce gradients of ice discharge, ice thickness, and surface velocity under debris cover. Independent of parameter selection, our simulations reproduce general relationships between debris cover percentages, AAR, and debris-perturbed surface velocity patterns from debris-covered glaciers.

The importance of debris flux and characteristic debris thickness on steady state glacier length
Increases in hillslope debris flux (ḋ flux ) lead to glacier extension (Figs. 8 and 9;Scherler et al., 2011b). But the rate and location of debris delivery to the surface will vary widely due to local geologic and climatic settings. Our simulations show that debris flux is more important in determining the steady state debris-covered glacier length thanḋ, d loc , or d width ( Fig. 9; Table 2). Processes of debris delivery to the glacier surface (e.g., deposition by avalanches, rockfall, the melt out of debris septa forming ice-stream interaction medial moraines, etc.) are first-order controls on the geometry of debris deposits on glaciers. Because debris flux trumps the importance ofḋ, d loc , and d width , the specific debris delivery pathway is secondary to the debris flux in determining glacier length at least for this 2-D case. The effects of changing h * are similar to the effects of varying the hillslope debris flux (Figs. 10 and 11). Establishing the importance of debris flux for individual glaciers requires that we constrain the variability of h * from glacier to glacier: small changes in h * can lead to large changes in steady state glacier length (Fig. 10). Simulations using an exponential debris thickness-melt curve (e.g., Konrad and Humphrey, 2000;Hagg et al., 2008) resulted in unrealistically long glaciers due to the rapid asymptote of melt towards zero (see Fig. 3). We argue that the hyperbolic parameterization (Eq. 3) is more physically defensible than the exponential, as we assume that heat transfer through debris is dominated by conduction.
Many paleoclimate estimates derived from glacial moraines neglect the potential effects of surface debris. Because debris strongly influences glacier length, independent of climate change, debris should be considered amongst temperature and precipitation as primary controls of paleoglacier lengths (e.g., Clark et al., 1994;Scherler et al., 2011b). The effect of debris on paleoclimate estimates can be minimized by avoiding de-glaciated catchments with high-relief headwalls, supraglacially sourced moraine sediments, or by using a debris-glacier-climate model to estimate the effect of debris on glacier extent.

The effect of steady debris input on patterns of Q, H and u surf
In all debris-perturbed simulations, the mass balance gradient down glacier from the location of initial debris emergence, x int , reverses relative to the debris-free profile, decreases toward zero, and becomes more uniform (excluding the terminal wedge; Fig. 5). This reversal results in a reduction of the surface mass balance b relative to the steady-state debris free glacier (Fig. 6). Reducing b toward zero reduces ice discharge gradients leading to glacier extension. Thick debris reduces b toward 0 and also makes b more uniform (Fig. 5). This leads to ice discharge gradients that are reduced toward zero and become more uniform near the terminus (Fig. 5). Because Q = H u, the surface velocity pattern follows a similar concave up pattern near the terminus where ice thicknesses are small and b is close to zero (Fig. 6). Low ice thicknesses and thick debris near the terminus leads to low, nearly uniform surface velocities, independent of climate change (Fig. 6). While it is possible that debris cover can produce low velocity portions of glaciers independent of climate change, periods of negative mass balance can also lead to extensive portions of debris-covered glaciers with low surface velocities due to the largest increases in melt rates occurring near x int (e.g., Kirkbride, 1993).
The ice discharge at the point of debris emergence, x int , controls the steady state glacier length and the down glacier patterns of ice discharge, ice thickness and surface velocity. In steady state, ice discharge at x int represents the volume of ice per unit time that must be ablated between x int and the terminus. Holding other debris-related variables constant, if debris emerges where ice discharge is high, the glacier will extend further because more glacier surface under thick debris (where melt rates are low and more uniform) is needed for ablation and match the large ice discharge at x int . If debris emerges where ice discharge is small the glacier does not extend as far because less area is needed under debris to match ice discharge at x int (Fig. 6). The location of debris deposition/emergence relative to the ELA is therefore an important variable in the debris-glacier system, as it controls the relationship between debris cover percentage, AAR, and the pattern of surface velocities (Fig. 11).
The specific terminal pattern of ice discharge and thickness is controlled by the rate of debris removal from the terminal wedge (Appendix A and B; Figs. A1 and B1). Ifḋ term flux is high, an ice cliff may persist at the toe leading to high melt rates and the pre-mature termination of a glacier when compared to a glacier with a lowḋ term flux . If the magnitude oḟ d term flux is low then the toe may be drowned in debris, and the glacier may never reach steady state even with a steady climate. The glacier would continue to accumulate debris and slowly advance down valley with a slightly positive net mass balance (e.g., Konrad and Humphrey, 2000). It may be useful to consider if individual debris-covered glaciers are accumulating debris mass through time, losing debris mass through time, or potentially in quasi-steady state with regard to debris (Fig. 4).
The response time of the modeled glaciers is therefore dependent on the parameterization ofḋ term flux (Appendix B). A glacier with rapid debris removal at the margin will tend to reach a steady state much faster than a glacier with slow debris removal from the margin (Appendix B). Documenting the rates of debris removal at the margin is therefore vital for modeling and understanding individual debris-covered glacier response.
In our steady state simulations, the ice thickness is increased up glacier from the point of debris emergence (Fig. 6). The thickness perturbations caused by emerging debris are diffused up glacier, leading to lower ice surface slopes and greater ice thicknesses than on debris-free glaciers forced by the same climate. The emergence of debris on a glacier can therefore perturb ice thickness both up and down glacier from the point of debris emergence. Debris cover decreases the surface mass balance and therefore also reduces the vertical component of englacial velocity; this leads to flow paths that are increasingly parallel to the surface (Konrad and Humphrey, 2000). Reducing ablation rates results in lower debris emergence rates, leading to the further advection of debris down glacier and expansion of the zone of debris emergence (Fig. 5a). Debris emergence zones will therefore tend to be wider than debris deposition zones.

Potential model improvements and future research
While we have explored first-order connections between glacier dynamics and debris deposition, additional components require investigation. Modeling the response of debriscovered glaciers to climate is the most pressing (e.g., Naito et al., 2000;Banerjee and Shankar, 2013;Rowan et al., 2015). The steady state results presented here can serve as initial conditions for future simulations exploring the response of debris-covered glaciers to climate change. Future efforts should further explore the importance of glacier size, environmental mass balance gradient, and valley bedrock profile as they modulate the effect of debris on glacier response.
We assumed a steady debris input for simplicity. In reality, hillslope erosion in high-relief settings occurs through thresholded, mass wasting processes. The effect of temporal and spatial changes in debris deposition must be addressed through both empirical and theoretical approaches. Isolated, large landslides have been shown to suppress melt rates, change glacier surface slopes and perturb glacier surface velocity fields (Gardner and Hewitt, 1990;Reznichenko et al., 2011;Shugar et al., 2012). If debris inputs are allowed to vary in space and time, a complex glacier length history will likely result even with a steady climate. The specifics of that history will depend strongly on the frequency and magnitude of mass wasting events and to a lesser degree the ice discharge at the point of debris emergence.
Our modeling did not account for the plan-view dimension of glaciers. Debris advected into the glacier between tributaries emerges to form ice-stream interaction medial moraines. While the spatial widening of such moraines has been addressed (Anderson, 2000), the merging of these medial moraines results in debris thickening that we do not account for. Our present work lays the framework for such a 2-D plan-view model.
Ice cliffs, surface ponds, and proglacial lakes are neglected in this study for simplicity but should be included in numerical models of glacier response to debris and climate change (e.g., Benn et al., 2012). Plan-view modeling of debris-covered glacier response is also needed (e.g., Menounos et al., 2013;Rowan et al., 2015). The melt-enhancing effects of thin debris covers should be included in future modeling efforts. Environmental mass balance profiles and snow lines are not steady from year-to-year. The response of debris-covered glaciers to interannual climate variability must also be explored (Roe and O'Neal, 2009;Anderson et al., 2014). Debris covers and glacier lengths will fluctuate in response to this variability due to the feedbacks between the debris emergence, ice dynamics, and climate.
Debris advection through and on a glacier can take hundreds of years, leading to memory in the system (i.e., the glacier responds to debris input from hundreds of years ago).
The response of individual debris-covered glaciers to climate change is therefore dependent on the distribution of debris on and in the glacier when the climate change occurs. Further constraint of englacial and surface debris is needed to predict the decadal to centennial response of present debris-covered glaciers to climate change.

Conclusions
It is necessary to constrain the effect of debris on glaciers so we can better predict the response of debris-covered glaciers to climate change. We provide a new framework to explore debris-covered glacier evolution and explore valley glacier sensitivity to debris input. Our simulations show the following: -For reasonable debris deposition fluxes, debris input can lead to glaciers that are many tens of percent longer than debris-free glaciers forced by the same climate but unperturbed by debris.
-Thick debris cover tends to reduce gradients of ice discharge, ice thickness, and surfaces velocities, independent of climate change.
-Debris-covered glacier length is highly sensitive to debris flux to the glacier surface. High surface debris fluxes can greatly increase glacier lengths relative to glaciers responding to the same climate without debris. Increases in debris flux lead to smaller AARs and larger debris covered fractions. Changes in the debris deposition zone width or the debris deposition rate are secondary to the total surface debris flux in governing the glacier geometry. This model provides a framework to quantify the effect of debris input on glacier length, and can therefore be used to estimate the effect of debris input on paleoclimate estimates derived from glacier models.
-The site of supraglacial debris deposition relative to the ELA modulates glacier response to debris. Steady debris input where ice discharge is high (near the ELA) leads to longer glaciers with greater fractional debris cover, whereas the same steady debris input where ice discharge is low (near the headwall or terminus) leads to shorter glaciers with smaller fractional debris cover.
-The importance of the mechanism of debris deposition onto glaciers (e.g., delivery by avalanching or by melt out of debris septa) is likely secondary to the importance of the total surface debris flux.
-Debris-covered glacier length is highly sensitive to the relationship between surface debris thickness and subdebris melt. Our simulations support the use of capped hyperbolic debris thickness-melt curve fits (Eq. 3) instead of exponential fits. -The rate and process of debris removal from the terminus exerts strong control on the time evolution of debriscovered glaciers, but only weakly influences the eventual steady-state length.
-Debris cover can perturb ice thicknesses and glacier surface slopes up glacier from the debris-covered portion of the glacier. Thick debris cover can expand the zone of debris emergence. Debris emergence zones will therefore be longer than zones of debris deposition.
Glacier response to debris cover is most sensitive to surface debris flux and the debris thickness-melt relationship. Our ability to predict the response of debris-covered glaciers to climate change, and to extract paleoclimate estimates from moraines in high-relief settings, is therefore highly dependent on our constraint of surface debris fluxes and debris thickness-melt relationship in the future and the past. Figure A1. The terminal wedge parameterization and debris removal from the model. Q in is the ice discharge into the terminal wedge.ḋ term flux is removed from the total volume of surface debris on the terminal wedge.

Appendix A
After the step change increase in debris deposition occurs, the steady-state debris free glacier evolves towards a debriscovered steady state. During this transition debris on the glacier surface is advected from cells with debris cover into debris-free cells. In our model, the debris thickness h debris (x, t) represents a layer of equal thickness on any cell. Debris thickens more slowly with a larger dx because the debris volume advected into a cell is spread over a larger area (due to the larger dx; dy=1; dy (in m)). There is therefore a timescale built into the thickening of debris in a cell that is dependent on dx. Because ablation rates are sensitive to debris-cover thickness, changing dx has an effect on glacier evolution. In order to test the effect of changing dx on the steady state debris-covered glacier length we increased dx from 100 (used in all simulations outside of this test) to 200 m. This test led to differences in steady state debriscovered glacier length which were less than 200 m even whenḋ flux was varied. The dx dependence does not effect the conclusions we draw from this study.
Without a terminus wedge parameterization, simulated glaciers advancing toward steady state become trapped in false steady states. Without a terminus wedge parameterization a new glacier cell is exposed to melt rates un-perturbed by debris. As a result, simulated glaciers become trapped in a steady length, even though large volumes of ice are melted without the protection of debris. To correct this, we implement a triangular terminal wedge parameterization for the last two grid points (the last ice-covered and the first ice-free grid point; Fig. A1; see Budd and Jenssen, 1975;Waddington, 1981) of the glacier which allows debris to cover the glacier terminus even when advancing or retreating. The volume and length of the terminal wedge is based on ice mass conservation. The volume of the terminal wedge at time t +dt is the sum of the old terminus volume, the ablated volume under debris, and the volumetric flow past the last grid point. Equation (16) and dx term , the surface length of the wedge, define the debris thickness on the terminal wedge.ḋ term flux removes debris from the total volume of debris on the terminal debris wedge. A single environmental melt rate is calculated based on the mean elevation of the terminal wedge, and subdebris ablation is calculated perpendicular to the surface of the wedge. When the terminal wedge length is greater than 2dx, the wedge parameterization moves to the next cell down valley. If the terminal wedge is shorter than dx the terminal wedge parameterization retreats one cell. Because the terminus parameterization allows the glacier to change length at the sub-dx scale, simulated glaciers avoid numerical traps and advance to true steady states. In this model, steady state occurs whenḋ flux =ḋ term flux and the glacier length is steady. Figure B1. Exploring various choices for theḋ term flux parameterization. Glacier lengths are normalized by the steady state debris free glacier length (L ssdf ). Irrespective of the choice of theḋ term flux parameterization the steady glacier length is nearly doubled. For all simulationsḋ flux is 3.2 m 3 m −1 yr −1 . The glacier will never reach steady state for choices whereḋ term flux cannot evolve to equalḋ flux . This occurs whenḋ term flux = c and c is less than 3.2 m 3 m −1 yr −1 . Circles represent simulations in which M surface (the total debris mass on the glacier) and glacier length did not reach steady state after 5000 years. The time labels show how long it took for the glacier to reach steady state for the cases whenḋ term flux = cḃh debris . All simulations presented outside of this plot use theḋ term flux = cḃh debris parameterization with c = 1 ( * in this figure).

Appendix B
Debris deposited on the glacier surface is removed from the glacier by ice cliff retreat or wasting down the terminal glacier slope. Unfortunately, the rates and processes of debris removal from glacier toes are poorly documented. We therefore explore parameterizations for the debris removal flux from the glacier (ḋ term flux ) and their effect on glacier length (using the base parameter set whereḋ flux = 3.2 m 3 m −1 yr −1 ). Each simulation starts with the ssdf glacier followed by a step change increase inḋ flux . We considerḋ term flux = c,ḋ term flux = ch debris , andḋ term flux = cḃ z h debris where c is a constant that ranges between 0.1 and 10 with variable units such thaṫ d term flux [in m 3 m −1 yr −1 ]. Independent of the parameterization, d term flux controls both the time needed to reach steady state as well as whether a simulated glacier can reach steady state (Fig. B1).
Large changes inḋ term flux lead to minor changes in glacier length even after 5000 years, implying that the choice of theḋ term flux parameterization would have a minor effect on the length results presented (Fig. B1). All three parameterizations lead to the same steady state length for low c values (190 % of L ssdf ).
Ifḋ term flux cannot evolve to a state whereḋ term flux =ḋ flux , surface debris thickens unrealistically and the glacier never reaches steady state. Forḋ term flux = c the glacier will never reach steady state if c is less than 3.2 m 3 m −1 yr −1 . Foṙ d term flux = ch debris , andḋ term flux = cḃ z h debris the value ofḋ term flux changes through each simulation based on the debris thickness on the toe and the local debris-free melt rate. Thė d term flux = cḃ z h debris parameter shows a wider length variation than theḋ term flux = ch debris parameterization becauseḋ term flux = cḃ z h debris results in a wider range ofḋ term flux values due to theḃ z term. To ensure that steady state can be achieved in each simulation, we include the melt rate term in theḋ term flux parameterization (Fig. B1) that codifies an assumption that debris removal processes at the toe are in some fashion dependent on local air temperature and hence melt rates. We useḋ term flux = cḃ z h debris for all simulations outside of this Appendix (with c = 1).