Numerical Simulation of Supraglacial Debris Mobility: Implications for Ablation and Landform Genesis

Supraglacial debris does not remain fixed atop ablating ice, but can move across the ice surface as supraglacial topography evolves. This active debris movement (distinct from passive movement due to underlying ice motion) affects landform genesis as well as the rate and spatial distribution of ablation. While observations of debris transport across evolving supraglacial topography are abundant, models of these coupled processes over timescales of decades and longer are few. Here I adapt a numerical model of coupled ablation and downslope debris transport to simulate the evolution of an idealized debris-covered glacier on the timescale of complete de-icing. The model includes ablation that depends on supraglacial debris thickness and a hillslope-scale debris transport function that scales non-linearly with slope angle. Ice thickness and debris distribution evolve with model time, allowing complete simulation of de-icing and landform construction in an idealized glacier test-section. The model produces supraglacial relief that leads to topographic inversions consistent with conceptual models of hummocky landform genesis. Model results indicate that the relief of the glacier surface and postglacial hummocks depend on the relationship between characteristic timescales for ablation and debris transport, which is defined as an index of debris mobility. When debris mobility is high, topographic inversions are rapid and supraglacial and postglacial relief are subdued. When debris mobility is low, more pronounced supraglacial relief is produced, but postglacial relief remains subdued. An intermediate mobility appears to optimize both postglacial relief and the rate of de-icingcompared with both highly-mobile and immobile debris. This enhancement of de-icing due to debris mobility could contribute to the observed anomalous rates of ablation in some debris-covered glaciers.


INTRODUCTION
Hummocky moraine landscapes are widespread along the margins of ancient ice sheets. A leading hypothesis for the origin of hummocky moraine invokes supraglacial debris movement during wastage of debris-covered ice (Sharp, 1949;Eyles, 1979;Clayton et al., 2008;Schomacker, 2008). According to this hypothesis, basins on the ice surface collect debris while inter-basin ridges shed debris. Once the ice is gone, thick debris from the supraglacial basins becomes hills, while areas of thin debris from surrounding ridgetops become swales and wetlands. This conversion of hills to basins during deglaciation is often called topographic inversion. At a 10 3 -10 5 m scale, hummocky moraine tracts often exhibit characteristic hummock wavelength and height (Johnson and Clayton, 2005). These characteristics must reflect the physical processes that shed or concentrate debris on an ablating, debris-covered ice surface.
Supraglacial debris not only serves as the medium for construction of distinct postglacial landforms, but also modifies ablation of the underlying glaciers. While thin (less than a few centimeters) or discontinuous debris cover can enhance melt rates compared with debris-free ice, thicker debris typically reduces ablation rates (Østrem, 1959;Nakawo and Young, 1981). The relationship between debris thickness and ablation rate is usually strong where ablation is dominated by downwasting, wherein heat for melting is transferred to the ice primarily by conduction through the debris layer. However, backwasting of ice cliffs, thermokarst-like melting from ponded and flowing supraglacial (and possibly subglacial) meltwater, and bottom melting can also contribute to ice mass loss to varying degrees (Driscoll, 1980;Schomacker, 2008;Thompson et al., 2016).
As recent climate warming has caused widespread glacier thinning and retreat, an increasing fraction of earth's glacier area is becoming debris-covered (Scherler et al., 2018), and this fraction is expected to grow in the near future (Herreid and Pellicciotti, 2020). While many studies have sought to simulate the physics of ablation on subannual timescales and local spatial scales in debris-covered glaciers, these models still perform inconsistently when scaled up to regional glacier mass balance assessments (e.g., Kääb et al., 2012). The difficulty may stem from inadequate accounting of unobserved processes such as changing extent of ponded (Thompson et al., 2016) or subsurface meltwater (Benn et al., 2017), ice dynamics (Banerjee, 2017), or elevation lapse rates (Vincent et al., 2016). However, debris mobilization can also influence the distribution of high-and lowablation-rate areas through time, a phenomenon that has yet to be explored in detail (cf. Nicholson et al., 2018). This paper explores some of the consequences of supraglacial debris mobility for the patterns and rates of ice ablation and debris distribution over multi-annual and longer timescales.
Though early descriptions of shallow slope failure and meltwater transport of supraglacial debris are abundant (Russell, 1893;Gilbert, 1904;Tarr and Martin, 1914;Sharp, 1949), not until the middle of the 20th century did researchers begin to explore the physics of debris movement (Sharp, 1949;Boulton, 1967). Although these authors discussed ideas about thresholds for debris movement and textural clues to transport mechanisms, formal mechanical descriptions of debris movement only emerged in the 1970s (Boulton and Paul, 1976;Lawson, 1979Lawson, , 1982Paul and Eyles, 1990). These authors recognized supraglacial debris flows of varying degrees of fluidity and sought explanations from frozen-ground engineering literature for failure thresholds and rates of movement. These discussions highlighted the role of meltwater pore pressure developed through "thaw consolidation" of supraglacial debris, which was reported to govern the rate and style of mass movement and the resulting sediment properties (Lawson, 1982) but was difficult to implement in practical models (Moore, 2018). Much more recently, Moore (2018) reframed the description of thresholds for debris destabilization in terms of material properties that are more readily quantified, and incorporated terrain attributes that allow prediction of the spatial distribution of unstable debris and zones of likely surface runoff. Observed relationships between supraglacial topography and debris thickness (Nicholson et al., 2018) and transient change in debris thickness across an evolving glacier surface  are broadly consistent with expected areas of debris destabilization. The stability model stops short, however, of predicting debris transport rates or deposition, and this remains an area of future research need.
Rapid mass movement and glaciofluvial transport are not the only mechanisms of downslope debris movement. Individual supraglacial debris particles have been observed to migrate relatively slowly across underlying ice in either the downslope or sunward direction (Fryxell, 1933). This typically occurs when the melt reduction beneath a particle causes it to become perched on a pedestal of ice relative to the surrounding ablating surface-features like this capped with large boulders are often referred to as a "boulder tables." As the pedestal height grows and/or the girth diminishes, the particle slides or topples from the pedestal to a new position where the process may begin again. The net result over repeated cycles is a relatively slow motion that some have termed "topple-walk." Fryxell (1933) found on Teton Glacier, United States, that most of the motion occurred in the downslope direction, though for gently sloping ice surfaces (< 12°) particles acquired trajectories increasingly deflected southward (toward the sun) with decreasing ice surface slopes. Anderson (2000) adopted a similar topple-walk mechanism in a model simulation of medial moraine evolution, rationalizing that downslope particle flux should be proportional to debris particle size and ice surface slope.
A recent study by Fyffe et al. (2020) explored topple-walk transport and related phenomena using high-resolution repeat imagery on Miage Glacier, Italy. The authors documented particle transport in areas of partial debris cover on the order of a few centimeters per day down slopes ranging from 5°to 30°. In the same areas, ablation proceeded at 3-6 cm per day, suggesting that downslope debris movement in this setting is of the same order of magnitude as ablation rate. This study largely confirmed the slope and particle-size dependence of debris flux by the topple-walk mechanism, but additionally documented slope-dependent creep in areas of complete debris cover at smaller (less than 1 cm/day) but still significant rates. Details of the creep mechanism were not investigated, but the process was an important component of the total debris flux across the ice surface.
In other terrestrial settings, creep of hillslope materials is caused by a variety of stochastic disturbances, both biotic and abiotic, which tend to slowly and intermittently displace particles downslope at rates that are sensitive to slope angle and slope curvature (Heimsath et al., 2002;Roering, 2004;Pawlik and Šamonil, 2018). Abiotic disturbances from these regimes that also act in supraglacial settings include wetting/drying cycles, freezing/thawing cycles and snow loading, among many others. Thus, while the physical processes responsible for topple-walk creep may not be sustained where supraglacial debris thickness exceeds a few clast diameters, several creep-inducing processes can be expected to be active. This remains an issue in need of further investigation.
In debris-covered glacier settings, it is important to realize that sub-debris ablation by downwasting requires debris to "settle" under the influence of gravity and perhaps an additional seepage force related to the evacuation of meltwater (Moore, 2018). Even if the circumstances of this seepage are not sufficient to surpass the stability thresholds described by Moore (2018), some downslope displacement may still occur. Much like the particle re-arrangements caused by acoustic vibrations introduced in some laboratory hillslope creep experiments (Roering et al., 2001), it is reasonable to hypothesize that the combined settling and seepage processes in supraglacial debris should promote net downslope transport that scales with slope angle (Furbish et al., 2009) and ablation rate (Houssais et al., 2021). This, too, requires future study. Nevertheless, it appears justified to proceed with a heuristic description of slow downslope debris transport until a more rigorous description, grounded in observation or theory, becomes available.

Coupled Melt-Transport Models
Evolution of debris distribution on debris-covered glaciers has been modeled in two different ways. At the scale of a full valley glacier, landsliding, debris emergence, melt-out, and passive supraglacial transport has been described by models seeking to explain. debris cover effects on glacier extent and moraine position (Shulmeister et al., 2009;Vacco et al., 2010;Anderson and Anderson, 2016;Anderson et al., 2018) or debris cover effects on mass balance and runoff (Rowan et al., 2015). In most of these models, debris is passively transported within or atop dynamic ice. A second type of model focuses on local-scale transport of debris across the ice surface as surface relief evolves through differential melting. This approach requires a description of active debris transport relative to the underlying ice, where any ice dynamics (and passive downglacier ice and debris transport) is treated separately. The local-scale coupling of ablation and debris transport in these types of models creates a complex pattern of evolving supraglacial relief that can govern both the spatial patterns of debris accumulation and the rate of ablation of the underlying ice. The present study builds on this latter type of model, so I review some key examples of this approach here. Anderson (2000) investigated how the width and topographic expression of medial moraines evolves as one moves downglacier from their englacial debris source near the glacier equilibrium line. Ablation under debris was approximated by an exponential decay function with increasing debris thickness, neglecting the small-thickness ablation enhancement characteristic of some Østrem curves. Spatial variation in debris cover thickness led with time to creation of relief on the glacier surface. Anderson then reasoned that debris was transported across the ablating ice surface by the topple-walk mechanism, with the probability distribution of toppling steps becoming larger and more downslope with increasing slope angle. He therefore reasoned that this process could be described more generally with a slopedependent flux of the form q s ∼ D∇z where q s is downslope sediment flux [L 2 T −1 ], ∇z is the debris surface slope, and D is a rate parameter. Model results simulating the ablation of a glacier cross-section with vertical englacial debris bands compared favorably with the patterns of medial moraine relief and width on an Alaskan valley glacier (Anderson, 2000).
Following in part from Anderson, Fowler and Mayer (2017) explored the origins of "ice sails," which are ridges or pinnacles of debris-free ice on the order of 10 0 − 10 1 meters high and 10 1 − 10 2 meters long that emerge in some settings from ice otherwise covered with thin debris. Their model focused on the physics of ablation that is enhanced (compared to debris-free ice) by the presence of a thin debris layer, which results in lowering of debris covered ice relative to clean ice. The model also included a description of diffusive (slope-dependent) debris flux similar to that used by Anderson (2000), which allowed emerging ice sails to shed debris to surrounding debris-covered ice lows. Their model produced results broadly consistent with observed ice sails on a Himalayan glacier (Fowler and Mayer, 2017). These features are, however, transient and unlikely to have significant impacts on either glacier mass balance or landform genesis. Nevertheless, they demonstrate one among the diverse array of features that can arise through coupling of ablation and local-scale debris transport.
A recent study by Mölg et al. (2020) explored the evolution of debris cover originating from medial moraines similar to those described in Anderson (2000). They constructed a model similar to that of Anderson (2000) except with a hyperbolic function of debris thickness describing sub-debris melt (Anderson and Anderson, 2016). The authors used this model to help explain both the evolution of supraglacial relief with downglacier position and the influence of that relief on supraglacial and subglacial meltwater drainage. They show that the results compare favorably with the topographic features of Zmuttgletscher in the Swiss Alps. They further link the development of ice cliffs, which contribute disproportionately to ablation in some debris-covered areas, to meltwater concentration in the troughs between moraine ridges.
The studies reviewed above demonstrate the power of simple models coupling ablation and debris redistribution to improve our understanding of surface features of debris-covered glaciers. Each of these studies, however, is focused on reproducing the features of a particular glacier (Anderson, 2000;Mölg et al., 2020) or generating a particular feature (Fowler and Mayer, 2017). None explores the range of features and patterns that emerge as model parameters vary, nor do they carry simulations to complete de-icing (de-icing is used here to refer to the long-term wastage of an ice volume of interest, independent of process or glacier dynamics). Thus, many process-form relationships remain unexplored, particularly on the longer timescales relevant to landform genesis.
In this paper, I combine some aspects of these prior studies to simulate the coupled processes of sub-debris ablation and debris transport in downwasting debris-covered glaciers. Attention is focused on the influence of initial conditions and parameter values on the coupled dynamics and consequent distribution of debris on the ice-free landscape at the end of simulations. Generality and simplicity is favored over detail and complexity in order to highlight the relationships that most strongly impact the morphological result.

MATERIALS AND METHODS
Ablation rate _ m under a debris layer is approximated here with a simple hyperbolic function of debris thickness, H: where K is thermal conductivity, ρ i and L f are the density of ice and latent heat of fusion for ice, respectively, and T s is the debris surface temperature. The first term on the right-hand side contains thermal properties of the system, while the second approximates an average linear temperature gradient through the debris layer. While temperature gradients within debris are rarely linear at any moment in time, the assumption of linearity is a convenient simplification that has also been substantiated in real-world settings (Nicholson and Benn, 2006;Rowan et al., 2021). The ablation enhancement at small debris thickness exhibited in some Østrem curves is omitted here for simplicity, as it has been in many other studies (Anderson, 2000;Anderson and Anderson, 2016). The parameter H p is a small reference debris thickness, which ensures that the hyperbolic function yields finite ablation rates in the limit of small H. This parameter also constrains the ablation rate for debris-free ice. When melt releases debris from debris-bearing ice, the resulting debris production P is described as where c is volumetric debris concentration in the ice and n is the porosity of the debris once it is released. Differential melt produces relief on the debris-covered ice surface and debris transport can take place by the various processes reviewed above. For the purposes of this analysis, transport is assumed to take place in one horizontal dimension x, governed by conservation of sediment mass: where as before H is the thickness of supraglacial debris, H z surf − z ice , and z surf and z ice are the elevation of the debris surface and debris-ice interface, respectively. Here, q x is the along-slope unit-width flux of sediment in the x-direction relative to the underlying ice, which does not move or deform. This conservation of mass rule states that the change in debris thickness at any place x on the ice surface is due to a difference in flux of sediment between x and x + δx plus any "production" of sediment P as debris is released from melting ice. In reality, the sediment flux q x depends on several processes that aren't easily summarized in a simple transport equation. However, on the long timescales captured by a landform evolution model, a heuristic transport equation that satisfies a few basic requirements may still yield illuminating results (Dietrich et al., 2003). A general form of nonlinear transport is used here to simulate the contributions of both slow creep and rapid mass movements to downslope debris flux. Transport rate thus increases sharply as the slope S zz surf /zx approaches a threshold, S c , above which shallow landsliding is presumed to occur: Here, D is a transport rate parameter analogous to diffusivity that constrains the volumetric debris flux per unit width for a reference slope. The parameter α affects the nonlinearity of the slope-dependence and is often taken to be equal to 2. Relationships of this form have been used extensively in modeling the degradation of mountainous, soil-covered landscapes (e.g., Roering et al., 2001). Anderson (2000) derived an expression for D that related downslope transport rate to the topple-walk magnitude and frequency. While the rationale for such a relationship may not hold for large debris thicknesses, a functionally-similar transport coefficient may be expressed: where D 0 is a reference rate constant (valid for large debris thickness) and H p remains a small reference debris thickness. When H is large compared with H p , D approaches D 0 , but for small H, D approaches zero. Like Anderson (2000), I make the simplifying assumption that gradients in debris thickness are negligible compared with slope gradients so that debris thickness change is dominated by slope curvature terms in the combined Eqs 3-5.

Modeling Approach
A 1D explicit finite difference model was constructed in MATLAB using Eqs 1-5 to explore the evolution of supraglacial topography and debris re-distribution, as well as the resulting patterns of debris accumulation following complete de-icing. A conservative finite-difference scheme adapted from Perron (2011) was used in which melt and debris production were computed at nodes, while slope and debris flux were computed at midpoints between nodes. The ablation and debris transport problems were each treated in one spatial dimension, where ablation caused adjustment of the ice surface vertically (z) and Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 710131 debris transport could occur across the ice surface horizontally (x). In each time step, ice and debris surface elevations were adjusted vertically according to computed ablation distributions, followed by computation of debris fluxes and thickness changes. While an adaptive time-step tied to the value of D would have ensured numerical stability, for simplicity and computational efficiency, time step length was fixed for each simulation to Δt aΔx/2D 0 , where Δx is the grid size (1 m in all simulations) and a is a manually-adjusted parameter < 1 that was tuned by trial-and-error to maintain numerical stability. Time steps ranged from approximately 0.01-0.5 years. Ice flow was not explicitly modeled, making simulations most relevant to stagnant ice. All simulations were run until all ice was gone from the domain (Code for this model will be freely available at https:// github.com/peteymoore/DCGsimulation through the author's Github site should this manuscript be accepted for publication).
The behavior of the model was investigated within two idealized domain types: 1) a rectangular glacier cross-section with one or more vertical englacial debris bands; and 2) a rectangular section with no englacial debris but a supraglacial debris blanket of non-uniform initial thickness (Figure 1). Noflux conditions were enforced on the left and right boundaries for both domain types, representing either barriers to transport there or a mirror-image (periodic) domain. Domain sizes and aspect ratios were varied along with parameter values to explore impacts on ablation and debris accumulation patterns. For simplicity, all parameter values were held constant during a given simulation, but varied among simulations. Parameter definitions and ranges of values used in simulations are given in Table 1.
Domain type 1 was first used to validate model results by comparison with results of Anderson (2000), and then to explore the effects of varying initial conditions (e.g., domain geometry and debris mobility) and physical properties on supraglacial relief and debris distribution during and after melt. In all domain 1 simulations, the englacial debris sources were arbitrarily defined to be 10 m wide. In most of these, a single debris band was located along the right-hand boundary of the domain as illustrated in Figure 1A. However, a few simulations were conducted with multiple englacial debris bands in the interior of the domain to explore the interactions between sources. Domain type 2,  representing debris originating from events like landslides that don't produce continuous debris sources with glacier depth, was used to further explore the effects of initial conditions and physical properties on long term melt rates and changes in debris thickness distribution. In all of the domain 2 simulations presented here, initial debris thickness ranged from 0.1 to 1.1 m, with an arithmetic mean initial thickness of 0.6 m. A smooth ramp in debris thickness spanning 10 nodes accommodated the transition between thick and thin debris ensuring that initial debris thickness gradients did not exceed 0.1. In each domain type, geometry, and parameter values were varied over reasonable ranges to explore consequences for deicing rates and surface relief. Ranges of parameter values considered are presented in Table 1.
The relationships between characteristic rates of transport and melt are expected to determine how supraglacial relief evolves with time, and therefore whether debris thickness evolves toward uniform thickness or growing non-uniformity. To explore this with more generality, the length dimensions can be scaled by a characteristic length of this system. With this in mind, results are compared across model runs using scaled spatial variables normalized to the arithmetic mean debris thickness H.
Assuming that melt occurs by conduction through supraglacial debris only, the timescale t m to melt a thickness ZH of ice under a uniform debris thickness H is ρ i L f ZH 2 /KT s . Here, H is taken to be the arithmetic mean thickness. For domain type 1, mean thickness is defined as the final debris area divided by domain width, while in domain type 2 it is the mean initial thickness. Since the evolution of debris thickness differs between these two settings, results are only compared within domain types. The timescale t t for debris transport over a distance XH is X 2 H 2 /D 0 , assuming that debris transport occurs primarily by creep (i.e., S < S c ). Define the mobility index M, as the ratio of t m to t t : Since each characteristic timescale is the reciprocal of a characteristic speed, M provides a means to characterize the relative speed of debris transport compared to supraglacial relief production by differential melt. For larger values of M, transport efficiently responds to any developing supraglacial relief, while for smaller values differential melt produces relief more rapidly compared with debris transport response. While this simple model cannot reasonably represent the full range of behaviors exhibited in supraglacial debris, it can nevertheless be used to explore the impacts of varying parameter values on surface relief and de-icing time relative to that assuming a constant static debris thickness. The results highlight the influence of debris mobility on differential ablation, surface relief, and landform genesis in debris covered glacier settings.

RESULTS
Simulations with medial moraines (domain type 1) were broadly consistent with prior model experiments (Anderson, 2000;Mölg et al., 2020). Figure 2 summarizes two simulations with different domain lengths X, thereby representing different transport timescales. Each line within a panel represents the ice or debris surface at a moment in simulation time (an isochron), with constant intervals between isochrons given in the caption. Animations of the first simulation as well as an example simulation from domain 2 are provided as Supplementary Material. In both simulations illustrated in Figure 2, debris melt-out supplies debris cover to the surface, producing surface relief due to subsequent differential melt. As ablation proceeds, the relief between the moraine crest and remaining debris-free ice grows until debris extends across the entire domain width, representing the merging of debris covers from adjacent moraines or deposition at a lateral glacier margin. The timing of this merging affects the subsequent rate of ablation and distribution of debris.
In Figure 2A, the debris sources are relatively close together (2XH), and the debris-covered moraine slope extends to cover the full domain width early in the simulation. As a consequence, the debris thickness in the basin between moraines grows. By contrast, the debris sources in Figure 2B are farther apart, and ice between the widening moraine ridges completely disappears before debris extends to cover the domain fully. As a consequence, limited debris thickening occurs in the basins between moraine ridges, and ablation proceeds more rapidly (Figure 3). This should be no surprise, as the same volume of debris is permitted to spread over a larger area of ice surface, attaining smaller mean thicknesses. Even so, it leads to the somewhat counterintuitive result that larger volumes (wider domains) of ice disappear faster than smaller when flanked by the same single debris band. For debris band spacing exceeding that in Figure 2B, there is no further change in ridge relief or ablation rates since debris-free ice beyond the ridge slope has disappeared before debris advances further.
As noted in Anderson (2000), the flanks of growing medial moraine ridges quickly evolve to a nearly constant slope just sufficient to deliver downslope the debris produced by ablation at the debris band, beneath the moraine crest. Steady moraine slopes ranged from 0.35 to 0.6 (approximately 20°-30°) for high to lowmobility cases, respectively. Simulations reported here are consistent with the importance of debris supply rate (ablation rate or debris concentration) and mobility in governing the steady-state moraine-flank slope. As supply rate increases or diffusivity decreases, the slope steepens. The effect of steepening moraine-flank slopes is similar to that of increasing moraine spacing, leading to narrower moraine ridges that permit debris-free de-icing to proceed longer.
Results of several simulations with domain type 1 suggest that the postglacial moraine relief, quantified as the interquartile range of debris thickness normalized to mean debris thickness (IQRnorm), is small in most cases. This is consistent with field observations that indicate poor preservation potential for linear englacial or supraglacial features like medial moraines (Evans, 2009). Even so, the results shown here indicate that relief is greatest in settings with intermediate mobility index M. Where M is large, debris accumulation in basins between moraines is rapid relative to ablation and differential melting leads to a reduction in moraine slope and reduced debris fluxes into these basins. Where M is small, on the other hand, debris redistribution occurs slowly compared to ablation and is distributed more uniformly at the ice-free toe of the moraine ridge as it decays.
The position of single, sustained englacial debris sources fixes the location of topographic highs, and consequently also the intervening lows where debris collects. However, when multiple bands interact, as illustrated in the three-band simulation in Figure 4 (with bands at 100, 200, and 300 m), greater debris accumulation in topographic lows receiving debris from multiple bands suppresses ice lowering there, flattening those areas in order to shed debris to remaining topographic lows. While each band initially produces a distinct moraine ridge (like the isochron highlighted in blue), the eventual merging of ridges and suppression of melt between them leads to a convex glacier cross-section (highlighted in red). The glacier surface becomes peaked at the middle band position and steepens beyond the adjacent bands in order to deliver their added debris flux toward the lateral margins.
Interesting patterns emerge in simulations in domain type 2 as well. Here, instead of a sustained englacial debris source, all debris in the system is prescribed as a non-uniform supraglacial debris blanket in the initial conditions, akin to a rockfall or landslide deposit. Figure 5 compares two simulations with the same initial conditions (debris thickness 1.1 m on the right half, 0.1 m on the left half) but different parameter values manifesting as different degrees of debris mobility. In the high-mobility case of panel (A), differential melt initially lowers the surface on the left-hand side of the domain where thin debris suppresses ablation less. Highly mobile debris quickly descends the resulting surface slope from right to left, accumulating on the left and suppressing ablation there. The slope reverses twice more before de-icing is complete, at which time debris is nearly uniformly distributed across the domain (IQRnorm 0.18). In Figure 5 panel (B), debris mobility is lower and greater relief (and a steeper slope) develops before significant transport takes place. Transport is sufficiently slow in this case that the left side of the domain is almost completely icefree before significant debris accumulation occurs there. Subsequent debris transport accumulates much of the debris to the left side of the domain, where it remains as a higherrelief (IQRnorm 1.19) postglacial mound following a single episode of topographic inversion. Figure 6 illustrates the effects of debris mobility on the ablation history for the simulation in Figure 5B. The heavy blue line shows the loss of ice cross-sectional area (or volume per unit distance normal to the study plane) as a function of simulation time. Also shown is the steady melt predicted for the same amount of debris cover if it were static and distributed FIGURE 2 | Glacier surface isochrones from simulations in domain type 1 with debris-band width 10 m, debris concentration c 0.1, T s 2°C, K 1 J/(s m°C), and D 0 1m 2 /yr. Each isochrone line in this figure (as well as Figures 4, 5) represents the debris-covered ice surface at an instant in time, separated by constant time intervals and advancing in time from top to bottom. The simulations are classified as low-and high-mobility based on differences in the domain width (A) High-mobility simulation (X/Z 1) with 5-years intervals between isochrones (B) Low-mobility simulation (X/Z 2) with 2.5-years intervals between isochrones.
FIGURE 3 | De-icing histories for the two simulations shown in Figure 2. The quantity of ice remaining is represented by the cross-sectional area of ice, since ablation rate and ice lowering are non-uniform across the domain. Solid blue line corresponds to de-icing from Figure 2A, and dashed red line to Figure 2B.
Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 710131 FIGURE 4 | Glacier surface isochrones at 25 years intervals from a 400-m wide, 100-m thick type 1 domain simulation with three debris bands centered at 100, 200, and 300 m distance. Band width was 10 m in each case. All other parameters were the same as the simulation shown in Figure 2. Blue and red isochrones highlight the differences between early surface evolution (blue) when each debris band has a prominent, separate medial moraine ridge, and later (red) when moraine ridges have merged and the ice surface has organized to shed debris to the lateral ice margins. Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 710131 8

Moore
Debris Mobility uniformly over the ice surface (black dash-dot pattern line). Finally, the dashed red line shows the expected melt history if the debris was arranged according to the initial conditions of Figure 5B but was unable to move laterally within the domain. The simulated duration of deicing is significantly shorter for mobile debris than for either of the other immobile-debris cases. Although the initial simulated area loss rate for mobile debris is intermediate between the immobile end-member cases, the mobile debris shedding from areas that experience periods of melt suppression ensures that ice throughout the domain is subject to periods of reduced melt suppression under thinned debris.
The effect of debris mobility on total de-icing time and postglacial moraine relief in ten domain type 2 simulations is shown in Figure 7. De-icing time (the time required to ablate all ice from the model domain) is normalized to the time required to remove all ice under a static, uniform debris layer of thickness H, the arithmetic mean thickness. Note that, as indicated in the immobile, non-uniform debris case (red dashed line) in Figure 6, the de-icing time for low-mobility debris is greater than that of the reference uniform-thickness debris. Similarly, extremely-high mobility debris tends to produce de-icing times that are asymptotic to the uniform-thickness reference, since highmobility tends to disperse debris more evenly across the ice surface. But an intermediate mobility leads to expedited deicing times, reaching almost 15% shorter than the reference case in the simulation shown in Figure 5B.
Postglacial moraine relief is also affected by debris mobility (orange symbols in Figure 7). Much like the effective de-icing rate, postglacial relief is greatest for debris of intermediate mobility. Highly-mobile debris tends to become distributed more evenly across the model domain, leading to smaller differences in debris thickness. Low-mobility (but still mobile) debris doesn't homogenize like highly-mobile debris, but produces somewhat lower relief. However, it should be noted that immobile non-uniform debris with the same initial distribution as in Figure 5 would retain an IQRnorm of 1.65, suggesting a more complex pattern of relief generation for vanishing mobility.

DISCUSSION
The model explored here greatly simplifies the physics of coupled melt and debris transport, and omits many phenomena known to be important in real debris-covered glacier settings. Even so, the simulations highlight relationships that likely contribute in real settings to ablation rates and debris accumulation patterns. A key FIGURE 6 | De-icing history for the domain type 2 simulation shown in Figure 5B (blue) compared with immobile debris of the same initial thickness distribution (red dashed) and uniform thickness H mean 0.6 m (black dash-dot).
FIGURE 7 | De-icing time and postglacial moraine relief as a function of debris mobility for ten simulations in type two domains spanning two orders of magnitude in D 0 . Simulated de-icing time normalized to de-icing time for uniform immobile debris is shown in blue, while moraine relief, quantified as the interquartile range (normalized to mean) of debris thickness following complete de-icing, is shown in orange. Debris mobility on the horizontal axis is expressed for convenience as the base-10 Logarithm of the mobility index M, given by Eq. 6. objective was to examine the effects of debris mobility on relief and debris accumulation patterns on downwasting ice. Simulation types representing debris release from englacial debris bands and supraglacial debris deposits allowed examination of ablation-debris transport feedbacks in two end-member scenarios common to today's debris-covered valley glaciers.
Results of simulations in type 1 domains are broadly consistent with observations of medial moraine emergence from englacial debris concentrations. Medial moraine ridges grow in relief and width as melt releases englacial debris that subsequently suppresses melt in the underlying ice. The relief thereby created promotes downslope movement of debris down the broadening medial moraine slopes by mechanisms that remain poorly understood, but that may include topple-walk mechanisms in thin debris covers. Indeed, sustaining expansion of medial moraine ridges requires a mechanism for slow transport of debris down slopes that are usually too gentle to produce rapid mass-wasting. The required transport rate should be governed by the rate of debris supply from upslope debris melt-out. While these phenomena remain poorly documented in real field settings,  described melt-season expansion of debris-covered moraine slopes over adjacent dirty ice areas where the debris wasn't evacuated by supraglacial streams.
Model results indicate that debris mobility, characterized here by M, affects the long-term de-icing rate and debris accumulation patterns in predictable ways. When manifested as differences in spacing of the debris-band sources, widely-spaced debris bands correspond (Eq. 6) to lower mobility than closely-spaced bands. Where bands are spaced sufficiently that much of the ice thickness in basins between moraine ridges is melted before moraine ridges merge, de-icing is efficient compared to cases where moraine ridges merge and produce thick debris accumulations between them. The optimal spacing appears to occur where the horizontal distance between debris bands normalized to ice thickness (X/Z) is just less than twice the reciprocal of the equilibrium moraine slope (tan θ, where θ is slope angle). The equilibrium slope is, in turn, a function of the debris supply rate and the transport rate parameter. Simulations with multiple interacting debris bands highlighted the influence of debris supply rate when the transport rate parameter was held fixed. When the glacier surface evolved to have a single central moraine crest (red isochrone in Figure 4), the glacier surface increased in slope at the position of the adjacent (100 and 300 m) debris bands, which added their contributions to debris supply, thus requiring a steeper moraine flank slope.
Domain type 2 simulations allowed additional exploration of the influence of initial conditions and debris mobility parameters on supraglacial relief and de-icing rate. High-mobility debris rapidly responds to relief production by differential melt, maintaining relatively small surface slopes, and low relief. The ice surface in these cases evolves through several reversals of topography and results in debris thickness distributions that approach uniform. Lower-mobility debris reduces the number of topographic reversals until only one reversal occurs, manifesting as complete topographic inversion.
An optimal geometry for rapid de-icing appears to exist for these simulations as well, and it may coincide with the conditions that produce a single complete topographic inversion. This arises when the debris mobility is just sufficient to allow nearly complete de-icing under thin debris before the wave of downslope-transported debris thickens debris and suppresses melt there. Not only is postglacial moraine relief greatest under these conditions, but de-icing rate is as much as 15% faster than it would be if the debris thickness was uniform. This relative acceleration of melt is presumably due to the redistribution of meltsuppressing debris to areas where ice has already been removed, as well as the nonlinearity of the sub-debris ablation function.
The extent to which debris mobility modifies time-integrated de-icing rates is likely underestimated in the present model for at least two reasons. First, retaining finite debris thicknesses is important for numerical stability in the model formulation used here. Therefore, debris thickness in the model remained at least 15-20 cm even on steeply-dipping slopes once they were covered with debris. In reality, areas of thinner debris (and therefore higher ablation rate) are under-represented in the model. A second reason to expect greater impacts from mobile debris than shown here is the neglect of the thin-debris enhancement of ablation that is sometimes documented. Since the sub-debris melt model employed here omits any such enhancement, allowance for this melt-enhancement effect when debris is thin would further elevate ablation rates on slopes steep enough to shed debris.
Another key element missing from the simulations presented here and in other published model studies is the removal of debris from the supraglacial setting by meltwater. Many recent observations highlight the importance of meltwater in maintaining high relief and steep slopes on debris-covered glaciers (e.g., Fyffe et al., 2020;Mölg et al., 2020). If meltwater removes mobile debris from the base of a debris-covered slope, accumulation and melt-suppression are inhibited there and the existing relief can be maintained or enhanced. This stands in contrast, however, to the phenomenon of topographic inversion, which is also widely documented in debris-covered glacier settings. Thus, further research is required to explore the conditions under which meltwater streams can allow topographic relief to persist on debris-covered glaciers, and where it may be unable to prevent topographic inversion.
Because the model simulations described here were idealized onedimensional experiments, the results don't readily compare with particular field settings. Surface topography on debris-covered glaciers and ice-free supraglacial landsystems is very much two-(or three-) dimensional, and a more sophisticated 2D model would be required to generate meaningful comparisons. Even so, some of the patterns observed here can help to inform interpretations of process-landform relationships. The length-scales used in the domain set-up (∼10 2 m) are generally consistent with those observed in topographic features of modern debris-covered valley glaciers (Bartlett et al., 2020), and for simulations that produced significant postglacial relief, moraine ridges are similarly spaced. This spacing is also consistent with some observed hummocky moraine tracts (albeit in only one dimension) (Johnson and Clayton, 2005), and is inherited from the prescribed spacing of debris source nonuniformity. Even so, simulations with high mobility or extremelylow mobility yielded very little postglacial topographic relief, which is consistent with past studies of controlled moraines (Evans, 2009).
Thus, an interesting consequence is that supraglacial debris redistribution can both create and destroy relief in de-glacial and postglacial settings, depending upon mobility.

CONCLUSION
Debris mobility on debris-covered glaciers can take the form of rapid mass movements, fluid-assisted flows, or slow creep. This movement is not only essential for genesis of key landforms in the so-called supraglacial landsystem (Johnson and Clayton, 2005), but affects the magnitude and spatial distribution of ablation. The model simulations presented here examine relationships between measures of debris mobility, subdebris ablation and the spatial distribution and redistribution of debris during and following de-icing. Among the key findings are that debris mobility can hasten ice mass loss compared to immobile debris, and that an optimum mobility exists where de-icing rate is greatest and postglacial moraine relief is high. This optimal mobility is related to the initial distribution of debris sources and the ice thickness to be melted. These results have implications for the practice of projecting mass balance in areas of debris cover, where predictions of mass loss often assume static debris or debris with minimal spatial variation in thickness. The results also inform interpretations of process-form relationships in supraglacial landsystems. Greater insights will, however, require a better understanding of the processes of slow downslope debris transport in supraglacial settings and the variables that govern them.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
PM is responsible for all aspects of this paper.

FUNDING
This study was funded by the U.S. National Science Foundation under grant EAR-1225880.