Modeling debris-covered glaciers : extension due to steady debris input

Introduction Conclusions References


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 Figures

Back Close
Full occurrence of debris-covered glaciers (e.g., the Himalaya and the Alaska Range; Scherler et al., 2011b).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 melt rate of underlying ice (e.g., Østrem, 1959;Shroder et al., 2000;Owen et al., 2003).If debris supply is high to a glacier surface, mass balance profiles can be greatly altered, leading to increases in glacier volume and length (e.g., Scherler et al., 2011b; Fig. 1).Thick debris cover on glaciers can also lead to low accumulation-area ratios (AARs; Scherler et al., 2011b).Paleoclimate estimates will be exaggerated if typical AARs are assumed when reconstructing past climate from former debris-covered glacial moraines.Debris-covered glaciers exhibit a wide range of responses to climate change (Scherler et al., 2011a).While Himalayan debris-free glaciers are almost coherently retreating, Himalayan debris-covered glaciers are not responding coherently to climate change.Some Himalayan debris-covered glaciers are advancing, others are stationary, and yet others are retreating (e.g., Raper and Braithwaite, 2006;Scherler et al., 2011a;Benn et al., 2012;Banerjee and Shankar, 2013).However, there is a strong trend toward negative mass balance for most of these debris-covered glaciers (Bolch et al., 2011;Benn et al., 2012).
In situ documentation of debris-covered glacier mass loss is made difficult by non-uniform debris thicknesses and the presence of scattered ice cliffs and surface ponds.As a result, in situ debris thickness, sub-debris melt rates, sub-debris ice thickness measurements and complete summer balances are sparse (WGMS, 2008).Measurements of englacial debris concentrations and distribution are yet more difficult to obtain, but vital for predicting debris-covered glacier response to climate change (e.g., Kirkbride and Deline, 2013).In addition, exploration of century-scale response of debris-covered glaciers to climate is limited by short satellite and observational periods (Bolch et al., 2011).Logistical realities therefore limit our ability to constrain Figures feedbacks between debris deposition rates, the englacial environment, the supraglacial environment, ice dynamics, and climate change.
While logistics limit our ability to directly observe some feedbacks, many of the most provocative conclusions relating debris and glacier response are based on remotelysensed data.Scherler et al., (2011b) provided an extensive inventory of remotelysensed velocity and debris coverage data from 287 glaciers in High Asia.They inferred that (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 debrisfree 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 also act as targets for debris-covered glacier models.
Numerical models can help quantify feedbacks within the climate-debris-glacier system (e.g., Konrad and Humphrey, 2000).Debris-covered glacier models have been used to explore the response of valley glaciers to (1) the constant 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).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 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 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).However, because of the complexity of the debris-glacier-climate system, many feedbacks remain unexplored.Introduction

Conclusions References
Tables Figures

Back Close
Full In this study, we isolate the effect of debris on valley glaciers independent of climate change.Debris fluxes, deposition rates, deposition zone widths, and deposition locations vary from glacier to glacier (Fig. 1), yet we know little about how changes in these debris related variables effect glaciers.So we ask: What about debris delivery controls glacier response?
In many debris-covered glacier systems, debris is deposited in the accumulation zone, advected through the glacier following englacial flowpaths, and emerges in the ablation zone (e.g., Boulton and Eyles, 1979;Owen and Derbyshire, 1989;Benn and Owen, 2002;Benn et al., 2012).In order to explore the response of glaciers to surface debris cover, we formulated a new transient 2-D numerical model (x, z) that couples debris deposition, englacial debris advection, debris emergence, surface debris advection, debris-melt coupling, and shallow-ice-approximation dynamics (Figs. 1 and  2).By coupling these components, we are able to explore the sensitivity of debriscovered glaciers to changes in debris input related variables (across the entire glacier) and compare our theory-based results with Scherler et al. (2011b)'s dataset.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.This study lays the foundation for future modeling efforts exploring the response of debris-covered glaciers to climate change.

Theory and numerical methods
We employ a 2-D finite difference numerical model (in downvalley and vertical, x and z) that can simulate the evolution of temperate valley glacier response to climate and debris.Forced numerical glacier models (e.g., Nye 1965;Budd and Jensen, 1975;Oerlemans, 1986;MacGregor et al., 2000;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, conservation of ice requires that where x is the distance along the glacier flowline, H is the local ice thickness, ḃ is the local specific balance, and This requires a prescribed mass balance field, and a prescription of the ice physics governing ice discharge.

Annual surface mass balance of ice in the absence of debris
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:

Conclusions References
Tables Figures

Back Close
Full 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.

Annual surface mass balance: effect of supraglacial debris
Sub-debris melt rate decreases in an exponential or hyperbolic fashion 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) and change based on 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, sub-debris melt rates approach bare-ice melt rates as debris thins (h debris h * ), and 6429 Figures

Back Close
Full asymptote towards zero melt as debris thickens (h debris h * ).We use h * values based on data from 15 studies (Fig. 3; h * = 0.066 ± 0.029 m (1σ), and ranges from 0.03 to 0.13 m).We also show the most likely exponential fit to the data for comparison to the most likely hyperbolic fit (Fig. 3).The exponential curve fit declines toward zero melt more rapidly than the hyperbolic fit.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 depth-averaged 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 [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 Introduction

Conclusions References
Tables Figures

Back Close
Full where 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, and provides a more conservative estimate of sliding velocities when τ b > τ c (Kessler et al., 2006).We have 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, 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 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.
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: employing the boundary condition that w = 0 at z = 0 (i.e., we assume no basal melt).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 either the upper glacier surface or at the glacier bed.Supraglacial debris deposition largely occurs by mass wasting from hillslopes above glaciers, while sub-glacial debris entrainment occurs through regelation and Figures

Back Close
Full 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 (European Alps); Owens and Derbyshire, 1989 (Karakoram); Ballantyne and Harris, 1994;Humlum, 2000 (Benn and Evans, 2010).For simplicity, we neglect englacial thrusting and icestream 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 (d flux [=] m 3 m −1 yr −1 ), the debris deposition rate ( ḋ ), the debris deposition zone width ( ḋwidth ), and the debris deposition location ( ḋloc ).In the model, d flux is representative of the integrated effects of ḋ and ḋ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 Figures 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, ˙ 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 viable deposition/hillslope erosion rates for high-relief mountain environments (e.g., Heimsath and McGlynn, 2008;Ouimet et al., 2009;Ward and Anderson, 2011).ḋwidth defines the 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 debris deposition location with the variable ḋloc , which we allow to vary from near the headwall to near the glacier terminus.ḋ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.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.
Once of change of concentration of debris within a parcel of ice is: where ).The first term on the right hand side represents the rate of change of C due to vertical strain of ice.Note that if the strain rate is negative, signifying vertical thinning of an ice column, debris concentration in the ice will increase.The second term represents the rate of change of C due to the longitudinal changes in glacier thickness.The third and fourth terms represent changes in C due to advection in the vertical and the horizontal directions, respectively.

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).The change of surface debris thickness with time at the glacier toe is:

Conclusions References
Tables Figures

Back Close
Full Varying this parameterization has a minor effect on glacier length, but can have a considerable effect on the temporal evolution of the glacier as d flux must equal d snout 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 the debris advection schemes.We also impose a two-step anti-diffusion correction algorithm to the advection scheme (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).The time step, dt, for ice-physics and debris advection is 0.01 years.
We select the base model parameters to represent the ablation zones of debriscovered glaciers in the Khumbu region of Nepal (Kayastha et al., 2000;Bolch et al., 2011;Benn et al., 2012;Shea et al., 2015).Base simulations run on a linear glacier bed with a basal slope of 8 % and a maximum bed elevation of 5200 m (Scherler, 2014).We use a d

Numerical experiments and results
We first demonstrate the transfer of debris between model components and demonstrate steady state.We then explore the differences between the ssdf glacier and debris-covered glaciers and explore relative importance of ḋ , ḋwidth , ḋloc , and d flux on glacier length.We then test the sensitivity of the model to changes in h * and φ.Last, we compare our results to data from real debris-covered glaciers in High Asia.

Demonstration of debris-covered glacier steady state and conservation of debris
In order to compare steady state glacier lengths between simulations with different d 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, 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 6437 Introduction

Conclusions References
Tables Figures

Back Close
Full 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 with a steady climate
We first highlight differences in length, Q, H, and u surf between the ssdf glacier and single simulated steady state debris-covered glacier (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).h debris increases down glacier from the point of initial debris emergence, ˙ x int , except near the glacier toe where the d snout flux parameterization reduces h debris (Fig. 5-6).Down glacier from ˙ 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 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).

Comparison of debris-covered glaciers with different debris input locations
Debris input location ( ḋ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.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).Another way of stating this: 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 d flux and ḋloc (Fig. 7).Changes in d flux are accomplished by changing ḋ with ḋwidth held constant.The importance of ḋloc on glacier length increases with larger d flux (Fig. 7).The pattern seen in Fig. 7 is insensitive to changes in the linear bed slope.Debris emergence/deposition at smaller Q leads to larger max(h debris ).Increasing d flux leads to increases in max(h debris ) and the percentage of the glacier covered with debris (Fig. 8).

Sensitivity of steady state glacier length to changes in debris deposition rate and debris deposition zone width
Increasing either ḋ or ḋwidth leads to increases in d flux , but their relative importance 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 ḋwidth on glacier length, we ran 36 simulations in which we vary ḋ , ḋwidth , and ḋloc .Steady state glacier length increases with ḋwidth when ḋloc and ḋ are held constant (Fig. 9c).
Steady state glacier length also increases with ḋ when ḋloc and ḋwidth are held constant Figures
The dependence of glacier length on ḋ and ḋwidth is not linear (Fig. 9).If we combine the effects of ḋ and ḋwidth by comparing the d flux with steady state glacier length we see that steady state glacier length is primarily dependent on d flux (Fig. 9e).Length enhancement by a factor of 2 or more is viable for the range of d flux explored.

Parameter sensitivity
We explore the sensitivity of the model to changes in h * and φ using the base parameter set for other parameters and inputs.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 from 140 to 250 % of L ssdf when h * is varied from the extremes of 0.0035 to 0.165 m (160-215 % for the 1σ range (0.037-0.095 m)).Glacier length is not as sensitive to the choice of debris porosity, φ (Fig. 10).Variation of φ between the extreme range of 0 and 0.45 leads to lengths that range from 160 to 195 % extension from L ssdf .

Comparison of model results with remote sensing derived data
Our model results show that steady, high debris fluxes onto glaciers lead to glacier lengthening and high percentages of debris cover (Figs.exposed bedrock in the headwall increase with steeper slopes, it follows that increased debris input onto the glacier should also increase glacier length and the percentage of the glacier covered with debris.Our steady state model results confirm this inference and show how changes in debris input variables can capture first-order trends from real debris-covered glaciers (Fig. 11).Scherler et al. (2011b) showed that large debris cover percentages correspond with small AARs outside the typical range of 0.5-0.7 seen on debris-free glaciers (e.g., Meier and Post, 1979).In our model simulations, increases in d flux lead to increases in both steady state glacier length, and fractional debris cover (Fig. 11a).With a fixed ELA, the AAR must therefore decrease with an increased d flux (Fig. 11a).Varying h * (using the base parameter set with no changes in d flux or ḋloc ; Fig. 10) has a similar effect to varying d flux (Fig. 11c and d).Changes in ḋloc lead to small changes in AAR but considerable changes in fractional debris cover (Fig. 11a).Scherler et al. (2011b) also showed that debris cover percentage correlated with the ratio of average u surf from the lower half of glaciers to the average u surf from the upper half of glaciers.Increasing d flux leads to lower u surf in the lower half of glaciers relative to u surf in the upper half of glaciers (Fig. 11b).Changing the location of debris input, ḋ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.
While the simulations plot within the data from Scherler et al. (2011b), our steady state model results do not account for the full data spread (Fig. 11a).Our ssdf glacier has an AAR of 0.5.Adding debris to the model only reduces AARs.Simulations with initial ssdf glaciers with higher AARs could reproduce more of the data.The Scherler dataset was collected from glaciers responding to periods of negative mass balance.Reduced surface velocities under debris cover (not necessarily stagnant) -resulting from debris-covered glacier response to climate change -could account for the data with low debris cover percentages and low ratios of half length mean ice surface velocities (Fig. 11b).Introduction

Conclusions References
Tables Figures

Back Close
Full Changing the linear bed slope leads to similar relationships between debris cover %, AAR, and surface velocity to the simulations using the base bed slope.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 ssdf 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).The specific relationship of glacier response to debris is therefore also dependent on glacier size, bed slope, and the environmental mass balance gradient.
This model-data comparison shows that viable changes in debris flux, debris deposition location, and h * can cause changes in debris cover percentage, AAR, and glacier surface velocities that correspond with patterns observed from real debriscovered glaciers.This lends support to the viability to our model framework, while also providing quantitative, theoretical support to previous data-based inferences.

Discussion
We explored the sensitivity of a new debris-covered glacier model to changes in various parameters and debris input related variables.Simulated glacier lengths are most sensitive to hillslope debris flux and the selection of the characteristic debris thickness.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 d snout flux , although steady state glacier length is not.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.Our model reproduces first-order relationships between debris cover percentages, AAR, and debris-perturbed surface velocity patterns from High Asian debris-covered glaciers.Introduction

Conclusions References
Tables Figures

Back Close
Full  et al., 2011b).But the rate and location of debris delivery to the surface ought to vary widely due to local geologic and climatic settings.Our simulations show that the flux of debris to the glacier surface, d flux , is more important in determining the steady state debris-covered glacier length than ḋ , ḋloc , or ḋwidth (Fig. 9).Debris delivery processes 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 d flux trumps the importance of ḋ , ḋloc , and ḋwidth , the specific debris delivery pathway may be secondary to the debris flux in determining glacier length.
The effects of changing h * are similar to the effects of varying the hillslope debris flux (Figs. 10 and 11).Establishing the importance of d 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 resulted in unrealistically long glaciers due to the rapid asymptote of melt towards zero.The hyperbolic parameterization is more physically defensible than the exponential parameterization if we assume that heat is transferred through debris by conduction.
Many paleoclimate estimates derived from glacial moraines neglect the potential effects of surface debris.Because debris can have a strong effect on 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 mitigated by avoiding de-glaciated catchments with high-relief headwalls and supraglacially sourced moraine sediments.Figures

Back Close
Full

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 toe cell; Fig. 5).This reversal results in a reduction of the surface mass balance b relative to the ssdf glacier (Fig. 6).Reducing b toward zero reduces ice discharge gradients.The glacier must extend in order to reach a steady state.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 = Hu, 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 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 et al., 1993).
The ice discharge at ˙ x int controls the steady state glacier length and the down glacier patterns of ice discharge, ice thickness and u surf .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 large, the glacier will extend further because more glacier surface under thick debris (where melt rates are more uniform) is needed to ablate 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

TCD Figures Back Close
Full 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 toe.If d snout 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 d snout flux .If the magnitude of d snout flux is low then the toe maybe 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 steady state with regard to debris (Fig. 4).
The response time of the modeled glaciers is therefore dependent on the parameterization of d snout flux (Appendix B).A glacier with rapid debris removal at the toe will tend to reach a steady state much faster than a glacier with slow debris removal from the toe (Appendix B).Documenting the rates of debris removal at the toe is 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, ˙ x int (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 of comparable sizes.The emergence of debris on a glacier can therefore perturb ice thickness both up and down glacier from ˙ x int .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 ice melt 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 on real glaciers will therefore tend to be wider than debris deposition zones.While we have explored first-order connections between glacier dynamics and debris deposition, additional components require investigation.Modeling the response of debris-covered glaciers to climate is the most pressing.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 also explore the importance of glacier size, environmental mass balance gradient, and bed slope as they modulate the effect of debris on glacier response.We assumed a steady debris input for simplicity.In reality, hillslope erosion in highrelief 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 planview 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 planview model.
Ice cliffs and surface ponds 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).Planview modeling of debris-covered glacier response is also needed (e.g., Menounos, et al., 2013).The melt-enhancing effects of thin debris covers should be included in future modeling efforts.Environmental mass balance profiles Introduction

Conclusions References
Tables Figures

Back Close
Full 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 because of 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.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 understand the decadal to centennial response of debris-covered glaciers to climate change.

Conclusions
Before modeling the response of debris-covered glaciers to a warming climate, it is necessary to constrain how debris effects glaciers -independent of climate change.We provide a new framework to explore debris-covered glacier evolution and explore valley glacier sensitivity to debris input.Our simulations show that: -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 Introduction

Conclusions References
Tables Figures

Back Close
Full 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 sub-debris melt.Our simulations support the use of capped hyperbolic debris thickness-melt curve fits (Eq.3).
-The rate and process of debris removal from the terminus exerts strong control on the time evolution of debris-covered glaciers, but only weakly controls 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 deposition zones will therefore be more narrow than zones of debris emergence.
Glacier response to debris cover is most sensitive to surface debris flux.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 in the future and the past.Introduction

Conclusions References
Tables Figures

Back Close
Full Appendix A In our model, the debris thickness h debris (x, t) represents a layer of equal thickness on any cell.There is therefore a timescale built into the thickening of debris in a cell that is dependent on dx.Increasing dx from 100 to 200 m leads to differences in steady state debris-covered glacier length that are less than 200 m even when debris flux is varied [=].Because melt (Fig. 3) is highly sensitive to debris thickness, a newly formed glacier cell at the toe can be exposed to melt rates un-perturbed by debris.As a result, the simulated glacier can be trapped in a steady length, although large amounts of ice are melted without the protection of debris.To correct this, we implement a triangular terminus parameterization (after Budd and Jenssen, 1975;Waddington, 1981).The volume of the terminal triangle at time t + dt is the sum of the old snout volume, the ablated volume, and the volumetric flow past the last grid point.A single environmental melt rate is calculated based on the mean elevation of the toe, and ablation is calculated perpendicular to the surface of the triangle.Equation ( 16) and the surface length of the wedge define the debris thickness on the snout.When the snout length is greater than  Full  (Østrem, 1959;Loomis, 1970;Khan, 1989;Mattson, et al., 1993;Lundstrom, 1993;Kayastha, et al., 2000;Lukas et al., 2005;Mihalcea, et al., 2006;Nicolson and Benn, 2006;Hagg, et al., 2008;Reid and Brock, 2010;Wang, 2011;Fyffe, 2012;Brook, et al., 2013;Anderson, 2014) (mean h * is 0.066 ± 0.029 m (1σ), and ranges from 0.03 to 0.13 m).These curve fits are used to determine the parameter ranges in Table 1 for h * .The best exponential fit is the mean of all the exponential curve fits; using sub-debris melt = ae Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 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-ice-approximation (SIA) and basal sliding parameterizations in Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | embedded in the glacier, C, the concentration of englacial debris [=] kg m −3 , will change only by straining of the ice.Taking an Eulerian point of view, the time rate Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | where d snout flux is the debris flux into the foreland from the toe [= ḃ 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).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 Introduction Discussion Paper | Discussion Paper | Discussion Paper | d flux = 3.2 m 3 m −1 yr −1 , ḋ = 8 mm yr −1 , ḋwidth of 400 m, ḋloc is 42 % from the headwall to the steady state length of the glacier, L ssdf .
Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | In order to highlight the effects of ḋloc on glacier length, Q, H, and u surf , we highlight three simulations where d flux = 3.2 m 3 m −1 yr −1 , ḋ = 8 mm yr −1 , and ḋwidth = 400 m are held constant between runs and ḋloc is varied.ḋloc is varied from near the top of the glacier (7 % from the headwall to L ssdf ; Figs. 5a and 6a and c), to near the ELA (42 % from the headwall to L ssdf ; Figs. 4, 5b, 6b and d), and near the debris-free glacier toe (98 % from the headwall to L ssdf ; Figs. 5c and 6c and e).
Discussion Paper | Discussion Paper | Discussion Paper | 8 and 9).Remote-sensing derived measurements of u surf and AAR provide insight into valley glacier response to debris.We compare our model results to Scherler et al., (2011b)'s inventory of 287 debris-covered glacier surface velocities, AARs, and debris cover percentages from High Asia.Scherler et al. (2011b) noted that debris cover percentage on glaciers correlates with steep above-glacier hillslopes.Because hillslope erosion rates and the percentage of Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | importance of debris flux and h * on steady state glacier length Increases in hillslope debris flux (d flux ) lead to glacier extension (Figs. 8 and 9; Scherler Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 2dx, the glacier advances to the next cell.If the snout is shorter than dx the glacier 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 d flux = d snout flux and the glacier length is steady.Appendix BDebris 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 (d snout flux ) and their effect on glacier length (using the base parameter set where d flux m 3 m −1 yr −1 ).Each simulation Discussion Paper | Discussion Paper | Discussion Paper | starts with the ssdf glacier followed by a step change increase in d flux .We consider d snout flux = c, d snout flux = ch debris , and d snout flux = c ḃz h debris where c is a constant that ranges between 0.1 and 10 and variable units such that d time needed to reach steady state as well as whether a simulated glacier can reach steady state (Fig. B1).Large changes in d snout flux lead to minor changes in glacier length even after 5000 years, implying that the choice of the d snout 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 d snout flux cannot evolve to a state where d snout flux = d flux , surface debris thickens unrealistically and the glacier never reaches steady state.For d snout flux = c the glacier will never reach steady state if c is less than 3.2 m 3 m −1 yr −1 .For d snout flux = ch debris , and d snout flux = c ḃz h debris the value of d snout flux changes through each simulation based on the debris thickness on the toe and the local debris-free melt rate.The d snout flux = c ḃz h debris parameter shows a wider length variation than the d snout flux = ch debris parameterization because d snout flux = c ḃz h debris results in a wider range of d snout flux values due to the ḃz term.To insure that steady state can be achieved in each simulation, we include the melt rate term in the d snout 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 d snout flux = c ḃz h debris for all simulations outside of this Appendix (with c = 1Discussion Paper | Discussion Paper | Discussion Paper |

Figure 2 .Figure 3 .
Figure 1.(a) Schematic of the debris-glacier system.Debris deposited on or emerging in the ablation zone reduces melt rates (above the critical 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.

Figure 4 .Figure 5 .Figure 6 .Figure 7 .Figure 8 .Figure 10 .Figure 11 .Figure 12 .
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 d flux = d snout flux and the glacier length is steady (see Appendix A).

Table 1 .
Parameters definitions and values.