Production and Transport of Supraglacial Debris: Insights From Cosmogenic 10Be and Numerical Modeling, Chhota Shigri Glacier, Indian Himalaya

Many mountain glaciers carry some amount of rocky debris on them, which modifies surface ablation rates. The debris is typically derived from erosion of the surrounding topography and its supraglacial extent is predominantly controlled by the relative accumulation rates of debris versus snow. Because Global Warming results in shrinking glaciers as well as thawing permafrost worldwide, changes in both rates will most likely affect the evolution of supraglacial debris cover and thus the response of glaciers to climate change. Here we report 10Be concentrations measured in five amalgamated debris samples collected from the main medial moraine of the Chhota Shigri Glacier, India. Results suggest headwall erosion rates that are ~0.5–1 mm year−1, and apparently increasing (10Be concentrations are decreasing) toward the present. We employed a numerical ice flow model that we combined with a new Lagrangian particle tracing routine to explore the impact of spatial and temporal variability in erosion rates and source areas on 10Be concentrations in the medial moraine. Our modeling results show that neither changes in source areas, related to the transient response of the glacier to ongoing climate change, nor four different scenarios of spatial and temporal variability in erosion rates could explain the observed trend in 10Be concentrations. Although not accounted for in our modeling explicitly, we suggest that the observed trend could be due to transiently enhanced erosion of recently deglaciated areas, or to greater spatial variability in erosion rates than explored in our models.


Introduction
Debris-covered glaciers are widespread within the Himalaya and other steep mountain ranges (Scherler et al., 2018). They testify to active erosion of ice-free bedrock hillslopes that tower above valley glaciers (headwalls), sometimes more than 1 km high. It is long known that thick debris cover significantly reduces surface ablation rates (Ogilvie, 1904;Østrem, 1959) and thereby influences glacial mass balances; but its dynamic evolution along with climatic and topographic changes is poorly understood. Observations from New Zealand (Kirkbride, 1993) and the European Alps (Deline, 2005) show that, following the Little Ice Age (LIA), some valley glaciers markedly increased their debris cover, thereby forcing a different and often slower response to warming compared to debris-free glaciers. Such increases in debris-covered areas are expected for debris-rich ice: ablation zones grow in response to warming and ice-flow rates decrease due to ice thinning (Banerjee & Shankar, 2013). Because most englacial debris is supplied from headwalls in the accumulation zones of these glaciers (Benn & Owen, 2002;Scherler et al., 2011), any change in the rate of debris supply could further affect the dynamic evolution of the debris cover, during both warming and cooling periods. Such changes might range from large individual rock avalanches that could trigger glacier advances (e.g., Tovar et al., 2008;Vacco et al., 2010), to changes in the frequency of rock falls. Because processes of headwall erosion likely depend on air temperatures (e.g., Gruber & Haeberli, 2007;Hales & Roering, 2007;Walder & Hallet, 1985), there may exist important links between climate change, headwall erosion, debris cover, and glacial response that have rarely been addressed so far (e.g., Scherler, 2014).
Better understanding of the coupling between ice-free bedrock hillslopes and glaciers in steep mountains requires means to assess headwall erosion rates. Some studies applied repeat laser scanning of steep rock walls to determine the amount of erosion (e.g., Abellàn et al., 2014;Hartmeyer et al., 2020). Although important, such studies typically cover small areas and short time spans, and may be biased by interannual variations and the effects of recent climate change. Furthermore, difficult access to many glacial headwalls limits the applicability of this approach. Recently, researchers have started to investigate the potential to estimate headwall erosion rates from the cosmogenic nuclide abundance in supraglacial debris (Heimsath & McGlynn, 2008;Seong et al., 2009). In a study on three valley glaciers in Alaska, Ward and Anderson (2011) measured cosmogenic-nuclide concentrations in samples collected from ablation-dominated medial moraines, i.e., longitudinal oriented bands of supraglacial debris that melt out of debris-rich septa that separate ice-flow units coming from different tributaries (Figure 1a). The benefit of this approach is that (1) the estimated erosion rates represent values that are spatially and temporally averaged across relatively large headwall areas; (2) the source areas of the debris can be reconstructed, e.g., from remote-sensing images; and (3) the sampling procedure is relatively straightforward. However, despite a promising initial application (Ward & Anderson, 2011), this approach has rarely been employed (e.g., Guillon et al., 2015;Sarr et al., 2019), and its applicability in different settings remains to be tested. Here we present headwall erosion rates derived from 10 Be concentrations in the ablation-dominated medial moraine of the Chhota Shigri Glacier, Indian Himalaya. We combine our empirical, field-based approach with a numerical model of headwall erosion and glacial debris transport of rocks with a new Lagrangian particle tracing routine. Our motivation for the numerical model is to test the sensitivity of debris 10 Be concentrations to changes in source areas that are related to the transient response of the glacier to climate change, and to four different scenarios of spatial and temporal variability in erosion rates.  10 Be concentration during transport of debris, due to exposure to cosmic rays. (c) Schematic illustration of temporal evolution of ablation-dominated medial moraine after model by Anderson (2000). Surface ablation progressively exhumes englacial debris that locally lowers surface melt rates and thereby creates relief, which induces lateral transport of debris.

Cosmogenic Nuclides in Ablation-Dominated Medial Moraines
Cosmogenic nuclides, such as in situ-produced 10 Be, are formed by the interaction of secondary cosmic rays with substrate near the Earth surface (e.g., Gosse & Phillips, 2001). These interactions attenuate the cosmic particle flux so that the production rate decreases approximately exponentially with depth into the ground. In erosional landscapes, rocks progressively approach the Earth surface at a rate that is set by the erosion rate of the surface (Lal, 1991). The 10 Be concentration of rocks at the surface is therefore a quantitative measure of the erosion rate, with low concentrations indicating high erosion rates and vice versa.
Converting measured 10 Be concentrations to erosion rates requires a model of the near-surface path of rocks or sediment during which 10 Be production occurred (Figures 1a and 1b). Rocks that erode from steep ice-free headwalls in the accumulation areas of glaciers fall down on the ice surface and are buried by snow, unless they fall directly into a bergschrund. The rocks thus become englacial debris. Ice flow then transports the debris into the ablation zone, where surficial melting reexposes it at the ice surface (e.g., Eyles & Rogerson, 1978). Much of the debris accumulates along the ice margin and is dragged along as lateral moraines. Two lateral moraines that merge at the confluence of neighboring ice units form a medial moraine. The debris exposed in medial moraines therefore constitutes an amalgamation of rock fall deposits from different source areas. This conceptual model of the pathway of debris through a valley glacier is supported by the widespread occurrence of medial moraines in glacial ablation zones worldwide.
Although there exist alternative pathways and source areas for supraglacial debris, these can be identified and excluded from sampling, based on simple observations. First, rocks that are eroded from the bed of a glacier can be brought to the surface by internal thrusting of the ice (e.g., Alley et al., 1997;Hambrey et al., 1999). Such processes are frequently observed on polar glaciers but less frequently on temperate valley glaciers. The resulting supraglacial deposit is aligned with the intersection of the thrust plane with the ice surface, which typically traverses the glacier perpendicular to the flow direction. It should also be noted that most theoretical models of subglacial erosion assume that erosion rates scale with basal sliding (e.g., Hallet, 1979;Iverson, 2012), which requires a slippery bed, achieved through basal melting, and which then leads to debris fluxes from the ice toward the bed, not the other way around. Second, rock falls may also originate from hillslopes or lateral moraines above the ablation zone of glaciers. However, these rarely travel far across the ice surface and are unlikely to mix with medial moraines (e.g., van Woerkom et al., 2019). In fact, the deposition of lateral moraines in the ablation zone often creates a topographic barrier that efficiently traps hillslope material in ice-marginal valleys. Third, large rock avalanches can catastrophically dump significant volumes of debris across glaciers in either the accumulation or ablation zones, but these are easily recognized in the field (e.g., Arsenault & Meigs, 2005;Shreve, 1966;Uhlmann et al., 2013).
The 10 Be concentration in rockfall debris exposed in medial moraines can be separated into (i) a part that was acquired during erosion of the bedrock-the part we are most interested in-and (ii) a part that was acquired during glacial transport. Using an idealized glacier model, Ward and Anderson (2011) have shown that production of cosmogenic nuclides during englacial transport is small for a wide range of glacial conditions and erosion rates on the order of 1-2 mm year −1 . There are two reasons for this: First, the rates of burial and reemergence from the ice are related to the rates of snow accumulation and ice ablation, which tend to be much higher than erosion rates. Second, most of the englacial transport is sufficiently deep within the ice where 10 Be production rates are negligible. Once the debris of the medial moraine is reexposed in the ablation zone, it is passively transported downglacier with the speed of the ice surface. In the absence of dynamic thickening or thinning due to velocity gradients, progressive ice melting results in thickening of the debris cover along the medial moraine ( Figure 1c). Because thicker debris reduces rates of ice melting under the debris (Østrem, 1959), medial moraines often develop a positive relief that induces lateral transport of the debris by downhill motion away from their center (Anderson, 2000). In that case, medial moraines get wider downglacier; debris that emerged higher up on the glacier moved outward and the core of a medial moraine exposes the most recently emerged debris ( Figure 1c). Finally, the additional accumulation of 10 Be in debris during its transport on the glacier surface is proportional to the surface production rate and the downglacier surface velocity.
Known challenges in this approach relate to the episodic nature of headwall erosion and the limited degree of mixing of sediments during glacial transport (Sarr et al., 2019;Ward & Anderson, 2011). First, steep bedrock walls in alpine environments typically erode by frequent and relatively small rock falls, but also by large and infrequent rock avalanches (e.g., Arsenault & Meigs, 2005;Matsuoka, 2008;Sanders et al., 2013).
Whereas cosmogenic nuclide concentrations in medial moraine debris can be related to erosion by rockfalls, it is complicated to do so for large rock avalanches (Ward & Anderson, 2011). Although large rock avalanche deposits on glaciers are easy to recognize in the field and can thus be avoided during sampling, they may nevertheless contribute significantly to eroding glacial landscapes. Neglecting these processes may thus lead to underestimating erosion rates. Second, mixing of sediment during glacial transport is limited, due to the laminar nature of ice flow. This fact can be seen as both a challenge and an opportunity. The challenge is that samples collected from small parts of medial moraines could be biased toward individual rockfall events. The opportunity is that samples collected along different medial moraines could provide a spatially more finely resolved picture of headwall erosion rates in glacial catchments. In addition, because glacial transport is slow and medial moraines are long, they offer the opportunity to obtain time series of erosion rates, provided that the gain and loss of cosmogenic nuclides during glacial transport can be quantified.
In summary, measurements of cosmogenic nuclides in medial moraines provide many opportunities to obtain valuable insights into the dynamics of glacial landscapes. However, the conversion of measured concentrations to erosion rates may not be as straightforward as commonly assumed for fluvial landscapes (e.g., Brown et al., 1995). In particular, limited mixing during glacial transport requires taking into account spatial variations in cosmogenic nuclide production rates, as these vary strongly with elevation. Moreover, the current shrinking of glaciers imposes a transience, which likely affects sediment transport trajectories and supraglacial debris cover. In order to move forward, yet honor the complexity of glacial landscapes, we have developed a numerical model that allows tracking of rock particles from their source areas through a glacier and concurrently recording the production and decay of cosmogenic nuclides. We have tested the utility of this model in conjunction with five samples we collected from the Chhota Shigri Glacier.

Study Area and Previous Work
The Chhota Shigri Glacier is located in the North Indian state of Himashal Pradesh, within a small basin that drains to the Chandra River ( Figure 2). The northerly oriented glacier has a length of~9 km, covers an area of 15.7 km 2 , and ranges from~4,200 to~6,100 m elevation (all elevation values in this paper refer to above mean sea level). Among several smaller tributaries it comprises two main branches that join at~4,800 m ( Figure 2a). The southeastern branch is approximately twice as large as the western branch. The resulting difference in ice discharge between these two branches shifts the medial moraine westward from the center of the trunk valley ice. The medial moraine originates at the base of a bedrock ridge that protrudes up tõ 450 m above the surrounding ice (Figures 3 and 2a), and which likely constitutes the main source area, constrained by tracing streamlines on the glacier surface (between subbasins 3 and 4 in Figure 4). The width of the medial moraine increases from~20 m near the junction of the two main glacier branches to~60 m at a distance of~3,000 m from the glacier terminus ( Figure 3). Closer to the terminus, the main medial moraine merges with several other but smaller medial moraines. In the field, we observed isolated boulders with diameters of up to several meters (Figure 2c), but most of the debris in the medial moraine is of centimeter to decimeter scale.
Glaciological mass balance measurements have been carried out since 2002 and indicate that the equilibrium line altitude (ELA), where annual mass loss equals mass gain, is currently at an elevation of 5,150 m (Wagnon et al., 2007). Specific annual mass balances have been largely negative with glacier-wide values of −1.4 to −1.2 m year −1 , although in the hydrological year 2004-2005, the mass balance was close to zero, with an ELA at~4,850 m (Wagnon et al., 2007). Recent mass balance modeling, based on downscaled global circulation model data as input, suggests that the ELA increased by more than 100 m between the periods 1955-1969and 2000-2014(Engelhardt et al., 2017. Precipitation records from nearby weather stations indicate that most snow accumulation occurs during the winter months (Azam et al., 2014;Wulf et al., 2010). Ground-penetrating radar (GPR) measurements were carried out in 2009 and allowed mapping of the ice-bedrock interface along five cross sections that are located between~4,400 and 5,000 m elevation (Azam et al., 2012). A complete map of estimated subglacial bedrock elevation was obtained by Frey et al. (2014) using the method developed by Huss and Farinotti (2012), and is used for the numerical modeling (see section 5.2). This method provides distributed estimates of ice thickness, based on the ice surface topography, an apparent glacier mass balance, and on ice flow mechanics. The bedrock in the catchment consists entirely of granitic rocks that belong to the High Himalayan Sequence (Steck, 2003).
Between the terminus of the Chhota Shigri Glacier and the Chandra Valley, prominent latero-frontal moraines document a more extensive glacier in the past (Ramanathan, 2011) ( Figure 2); but we are unaware of any age dating of these deposits. However, surface exposure ages from erratic boulders and glacially polished surfaces in the Chandra Valley document an extensive valley glacier that retreated past the confluence with the Chhota Shigri tributary, about 17,000 years ago (Eugster et al., 2016). Other dated glacial deposits in the Chandra Valley provide evidence that many glaciers were more extended than today during the early Holocene, and at least in one case~3,000 years ago (Owen et al., 2001). It is therefore difficult to estimate View is to the south. (c) Steep bedrock headwall near the southern-most limit of the glacier. (d) Medial moraine and debris-rich ice exposed at a crevasse that is oriented transverse to the ice flow direction. View is to the north. the actual age of the latero-frontal moraines in the Chhota Shigri basin and to precisely resolve the glacial chronology. However, based on the reconstructed LIA extent of the nearby (5 km to the east) Bara Shigri Glacier (Chand et al., 2017), and the widespread occurrence of LIA-moraines throughout the Himalaya (Rowan, 2017), we suggest that these moraines likely correspond to the LIA.

Sample Collection and Preparation
We collected five samples from the main medial moraine of the Chhota Shigri Glacier between 4,638 and 4,773 m elevation. The sample locations are separated by approximately 500 m along the glacier ( Figure 3). All samples consist of an amalgamation of >1,000 surface clasts with grain sizes betweeñ 1 and~30 mm, yielding samples of~6 kg weight. The clasts were randomly taken from surface patches of approximately 5-10 by 5-10 m, which were centered on the medial moraine. All samples were processed in the Helmholtz Laboratory for the Geochemistry of the Earth Surface (HELGES) at the German Research Centre for Geosciences in Potsdam, Germany, following procedures described in von Blanckenburg et al. (2004), and briefly summarized here. The samples were crushed and sieved and we retained the 500-1,000 μm-sized fraction for further analysis. After standard physical separation of quartz, carbonates were dissolved in hydrochloric acid, followed by repeated etching with a mixture of weak hydrofluoric and nitric acids in heated ultrasonic baths. Quartz purity check was done on~2-g aliquots with ICP-OES. The samples were spiked with a 372-ppm 9 Be carrier, dissolved in concentrated hydrofluoric acid and subsequently evaporated. Newly formed silica-fluorides were broken up with aqua regia. Beryllium was extracted by ion exchange chromatography, precipitated from solution at pH 10, and dehydrated at 1000°C before packing into targets for accelerator mass spectrometry measurements at the University of Cologne, Germany (Dewald et al., 2013). 10 Be/ 9 Be ratios were determined relative to standard 07KNSTD with a nominal 10 Be/ 9 Be ratio of 2.79 × 10 −11 . A process blank with a 10 Be/ 9 Be ratio of 9.41 × 10 −16 was subtracted from all samples. Oblique aerial view of the medial moraine of the Chhota Shigri Glacier from Google Earth (https://earth. google.com). View is to the south and upglacier. Yellow pins indicate sample positions (from CSM1 in the distance to CSM5 in the front). Note spatially restricted patches of debris cover near uppermost samples, presumably related to landslide deposits.

Cosmogenic Nuclide Production
In isotopic steady state, i.e., when the in situ-production of cosmogenic 10 Be equals its removal by decay and erosion, the 10 Be concentration N at depth z can be described by (Lal, 1991) where P 0 is the surface production rate by neutron spallation, λ is the decay constant (based on a half-life of 1.36 Ma; Chmeleff et al., 2010;Korschinek et al., 2010), E is the erosion rate, ρ is the bedrock density, and Λ is the absorption mean free path. Production rates were calculated with the spallation production rate-scaling model of Lal (1991), as modified by Stone (2000), and a reference sea-level high-latitude production rate of 4.01 atoms (g year) −1 (Borchers et al., 2016). We did not consider muogenic production, because it is small (<2%) compared to spallationinduced production, whereas the geomorphic uncertainties are significantly larger (cf., Ward & Anderson, 2011). We used the numerical functions of the CRONUS online calculator v2.3 (Balco et al., 2008) to compute a grid of production rates based on an SRTM 1-arc second (~30 m resolution) digital elevation model (DEM). We used the formulas given in Dunne et al. (1999) and the TopoToolbox v2 (Schwanghart & Scherler, 2014) to compute shielding of cosmic rays by the surrounding topography. We also corrected the mass attenuation length as a function of the surface slope following Masarik et al. (2000).

Sample Results
Our five samples yielded 10 Be concentrations of 27-62 × 10 3 atoms g −1 ( Table 1). The 10 Be concentration of the upper four samples increases systematically in the downglacier direction, whereas the lowest sample does not follow this trend ( Figure 5). Assuming that the sampled debris originated from the headwalls that constitute the southern boundary of the Chhota Shigri Glacier (Figure 4), we obtain slope-perpendicular, ice-free headwall erosion rates of 0.6-1.3 mm year −1 (Table 1).

Potential Reasons for Downglacier Changes in 10 Be Concentration
An interesting question that arises from our sample results is what mechanisms may account for the observed downglacier increase in the 10 Be concentration? Because our samples all consist of thousands of individual grains of similar grain sizes, we do not think that the observed trend results from a bias due to unequal weighting of certain grains in the amalgamated samples. However, we acknowledge that the limited number of samples, together with the fact that only four of the five samples follow this trend, may raise the possibility that it is an artifact, borne out of random scatter. We do not want to disregard this possibility, but there also exist a number of processes that could potentially account for such a trend, and we address these in the following.

Accumulation During Surface Transport
Debris that has emerged from depth is passively transported downglacier at the speed of the ice surface. During transport, the debris further accumulates 10 Be, which may thus account for a downglacier increase in concentration ( Figure 1). The surface production rate of 10 Be at the mean sample elevation (4,715 m) is 67 atoms (g year) −1 . To account for the entire observed increase of~35,000 atoms g −1 over a downglacier distance of 1,500 m, the ice surface velocity would have to be~2.5 m year −1 . This is in contrast to an observed surface velocity of~30 m year −1 (Azam et al., 2012;Scherler et al., 2011;Wagnon et al., 2007), which would only account for an increase of~3,300 atoms g −1 (<10%) over the same distance. The actual 10 Be accumulation in our samples during surface transport is most likely even lower than the above estimate for two reasons. First, in an ablation-dominated medial moraine that is continuously fed with englacial debris, the supraglacial debris is transported away from the center of the medial moraine and is replaced by new englacial debris that emerges farther downglacier (Anderson, 2000) ( Figure 1c). Therefore, we cannot simply substitute space for time and assume we sampled similar deposits on their way toward the glacier terminus. Instead, each of our samples likely comprises a mixture of debris that emerged at the surface more or less recently, and which followed different trajectories during englacial transport. Second, the Chhota Shigri Glacier has been shrinking and thinning over the last couple of decades (e.g., Engelhardt et al., 2017), probably from a more extensive size during the LIA. Because the ice surface velocity generally increases with ice thickness, the glacier has most likely been faster than now in most of our samples' transport period.

Different Source Areas
The measured 10 Be concentration is a function of the 10 Be production rate and the erosion rate in the source areas of the debris. Because production rates increase with elevation, the observed increase in the 10 Be concentration could theoretically simply reflect higher-lying source areas of samples collected farther downglacier. 10 Be production rates in the source areas range between 54 and 93 atoms (g year) −1 , but 90% of the production rates range between 63 and 84 atoms (g year) −1 . If erosion rates in the source areas were similar, the latter difference in production rate would account for no more than 16,000 atoms g −1 difference in concentration (~45% of the observed increase). However, this number should be considered a maximum value that is unlikely to be reflected in our samples for two reasons. First, each sample consists of an amalgamation of many grains that likely stem from different source areas. This is particularly relevant for medial moraines, as they comprise two merged lateral moraines ( Figure 1a). Second, source areas are typically steep headwalls that extend over considerable elevation ranges (e.g., Figure 3). Therefore, even spatially distinct locations at the edges of the glacier surface are supplied with debris from different elevations and thus production rates.
A complicating factor is the prevailing imbalance of the glacier (Azam et al., 2014), most likely since the LIA (e.g., Chand et al., 2017), which accounts for glacier shrinking and associated changes in debris trajectories and available source areas. The source areas of our samples may therefore have changed through time due to the glacier's response to past climatic changes. For example, our lowermost sample could be composed of rocks that were eroded during the LIA from headwalls that are different from those contributing to the medial moraine at present (Figure 3).

Changes in Headwall Erosion Rate
A remaining possibility are spatial and/or temporal changes in headwall erosion rates. As noted earlier, Anderson's (2000) model of ablation-dominated medial moraines suggests that they expose debris, which has been delivered to the glacier at different positions in the accumulation zone (cf., Ward & Anderson, 2011). Debris that emerges first is expected to have the shortest travel distance and thus englacial transport times, whereas debris that emerges close to the terminus would have longer travel distance and transport time. For these reasons, the gradual increase in 10 Be concentrations could reflect (i) increasingly lower erosion rates higher up in the catchment, or (ii) increasingly higher erosion rates Headwall retreat rates were calculated based on an average source area production rate P 0 = 72.56 at/(gqz year) −1 , and an average attenuation length L = 137.67 g cm −2 .

10.1029/2020JF005586
Journal of Geophysical Research: Earth Surface toward the present. The likelihood of these explanations depends on (i) how well we can link our samples to individual source areas and how different the source areas are, e.g., how steep they are, and (ii) how different the englacial transport times are for our samples, and if these times are sufficient for changes in erosion rates to be recorded by 10 Be concentrations (e.g., Schaller & Ehlers, 2006). The possibility of spatial and/or temporal changes in headwall erosion rates may also be combined with temporal changes in source areas, due to the glacier's response to past climatic changes (see section 4.4.2).

Summary of Sample Results
The above discussion shows that several processes exist, which can explain the observed increase in 10 Be concentrations, but that it is difficult to evaluate them. Furthermore, despite the apparently simple glacier and catchment geometry, the time-transgressive nature of the exposed supraglacial debris requires an assessment of how past and present glacier changes play out in the production and transport of debris within the catchment. For the remainder of the paper, we seek to explore how sensitive the 10 Be concentrations along the medial moraine of the Chhota Shigri Glacier are to these factors and whether the observed trend could possibly be related to one or several of them. We do this with the help of a numerical model that explicitly accounts for the transport of rock particles through the ice and for the in situ-production and decay of cosmogenic 10 Be. We apply the model to the Chhota Shigri Glacier here, but it can in principle be applied to other glaciers, too.

Model Description
We used the coupled ice and landscape evolution model iSOSIA (Egholm et al., 2011), and implemented a new Lagrangian particle tracing routine to model the evolution of 10 Be concentrations in glacial debris. With iSOSIA, the flow of ice across an arbitrary land surface is computed using a second-order shallow ice approximation (Hutter, 1983) that includes higher-order terms related to longitudinal and transverse stress components. A detailed description of the model is provided in Egholm et al. (2011Egholm et al. ( , 2012. Here we focus on the new model component that allows transport of particles by the ice, while tracking their cosmogenic nuclide concentration . The particle tracing has two main advantages over a Eulerian implementation of glacial debris transport (cf., Rowan et al., 2015). First, numerical diffusion during the passive transport of particles by the ice is reduced. This point is particularly important when studying discrete features like medial moraines ( Figure 6). Second, each particle can be assigned properties, such as cosmogenic nuclide concentrations, that are transported with the particle and possibly modified during the simulation.
A sediment particle is formed when the sediment thickness in a grid cell exceeds 1 mm. The particle then collects the sediment and it starts to move. As the size of grid cells is 25 × 25 m, every particle carries a minimum of 0.625 m 3 of sediment. The exact starting position within the grid cell is chosen randomly from a uniform likelihood distribution. In ice-free parts of the landscape, the particle velocity u ! p is assumed to align with the steepest descent direction: where K h is a constant and h is the bed elevation. This is a simple linear diffusion model. Its purpose is to move particles formed on hillslopes and headwalls down to the ice surface below.
If a particle is formed in, or later moves into, an ice-covered cell, the particle velocity is computed by interpolating the local ice velocity to the current position of the particle. If a particle moves from an ice-free cell into an ice-covered cell, we assume that the particle starts at the ice surface. Instantaneous burial to greater depth by means of crevasses or bergschrund is not considered. The ice velocity is a three-dimensional vector, including vertical components arising from surface accumulation/ablation, basal melting, and horizontal compression/extension. iSOSIA is a higher-order shallow-ice model solving for depth-averaged velocities. However, in order to include velocity variations at depth within the ice, we adopted the standard assumption for shallow-ice approximations that the horizontal ice velocity caused by viscous ice deformation decays as a fourth-order polynomial down through the ice. The particle velocity parallel to the bed is then where u ! d is the depth-averaged ice-deformation velocity vector and u ! b is basal sliding; H is ice thickness and z is the burial depth in the ice (z = 0 at the ice surface). The velocity vector is transferred from the staggered-grid ice-velocity nodes to the position of a particle using bilinear interpolation. The burial depth in ice changes in time due to the ice mass balance: where _ M s is surface accumulation (positive) or ablation (negative), _ M b is basal melting (positive), and _ ε zz is the vertical strain rate (positive for ice thickening and negative for thinning). The vertical strain is generally a consequence of horizontal extension and compression. The particle position and burial in the ice is updated using forward Euler time-integration. Particles that are at the ice surface influence ablation rates according to (Anderson & Anderson, 2016) where _ M s0 is the clean-ice melt rate, h debris is the debris thickness, and h * is a characteristic length scale of the melt-damping effect of debris cover (Østrem, 1959) that we set to 0.1 m. The average debris thickness per grid cell is computed from the combined mass of all particles at the surface. Supraglacial particles were identified as particles with a burial depth of less than 0.1% of the ice thickness.
When a particle is formed from bedrock erosion it is assigned an initial nuclide concentration computed from Equation 1, assuming steady state balance between bedrock nuclide concentration and local erosion rate. The nuclide concentration of particles is updated through time using where N t and N t+dt are nuclide concentrations before and after time step dt and P(z) is the nuclide production rate at burial depth z, computed from where ρ i = 900 kg m −3 is the ice density.

Model Setup and Boundary Conditions
We used the estimated ice thickness maps from Frey et al. (2014), subtracted from our 30-m resolution DEM, to obtain the bed topography beneath the present-day ice cover (see sections 3 and 4.2; Figure 7). The ice thickness maps are based on the approach of Huss and Farinotti (2012), and glacier outlines from the Randolph Glacier Inventory (RGI Consortium, 2017). However, the outlines of the Chhota Shigri Glacier have been modified in between different versions of the inventory, and some of the tributary glaciers are not connected to the main trunk glacier of Chhota Shigri (c.f. Figure 4c in Huss & Farinotti, 2012). We used high-resolution images in Google Earth (https://earth.google.com/web/) to assess actual connections between different ice bodies and manually merged the ice thickness maps based on version 3.2 and 5.0 of the RGI. The estimated ice thickness compares reasonably well with GPR profiles (Azam et al., 2012), except for the profile close to the terminus, where the mismatch is greater (Figure 7). However, as our samples are located upglacier from this area, it does not affect our results and conclusions.
We simulated ice ablation using a simple positive-degree-day approach, in which we tuned the parameters to get an approximate match between the modeled and observed vertical mass-balance profiles (Wagnon et al., 2007). For the positive-degree-day model we assumed the atmospheric lapse rate to be 6°C km −1 and the seasonal temperature variation to be 10°C (e.g., Scherler, 2014). The maximum amount of ice accumulation was also guided by the mass balance measurements, whereas the pattern of accumulation was adjusted based on the modeled ice discharge from the eastern and western branch of the glacier and its tributaries. In particular, we included a north-south decrease in ice accumulation, which follows the strong precipitation gradient in the study area (e.g., Wulf et al., 2010). We also included snow avalanching to transport snow from steep hillslopes (slopes > 25%) above the glacier to the glacier surface (Kessler et al., 2006).
The spatial domain of the model comprises the catchment of the Chhota Shigri Glacier, which we discretized into a staggered grid of 25-m cell size. iSOSIA uses adaptative time stepping based on the Courant-Friedrichs-Lewy criteria for numerical stability (Δt ≤ Δx 2u max , where Δx is cell size and u max is the maximum ice velocity in the grid) resulting in time steps varying between 0.01 and 0.1 year. After a spin-up time of 1,000 years, when ice volume reached a steady state, we ran our simulations for a time period of 3,000 years, in which we varied only atmospheric temperatures (Figure 8). We started the simulations with our climate parameters tuned so that the initial steady-state ice extent mimics the present-day ice-extent. A −0.6°C step decrease in atmospheric temperature at 2,000 years induced glacier advance to a position close to the latero-frontal moraines that we relate to the maximum LIA extent (Figure 2a). The glacier advance was reversed (ΔT = +0.6°C, relative to the initial increase) at 2,500 years, and the glacier was forced to further retreat at 3,000 years, when the atmospheric temperature was increased once more by +0.7°C. We note that the above climate scenario is not meant to match reality, or to provide a prediction of expected future changes. It rather serves the purpose of triggering glacial changes and to examine their impact on 10 Be concentrations in supraglacial debris. We selected four points in time that we will repeatedly come back to in the following. These times coincide with the three step changes in temperature (at 2,000, 2,500, and 3,000 years), and a fourth at 3,500 years, when the glacier has largely adjusted to the last temperature change and achieved its smallest extent during the simulation.

10.1029/2020JF005586
Journal of Geophysical Research: Earth Surface

Sediment Production
To explore the above-mentioned potential reasons for the observed downglacier increase in 10 Be concentration, we conducted four numerical experiments that differed only in the way sediment is produced at ice-free headwalls (Figure 9). In the first model (Uniform), we set the rate of sediment production, E, to a uniform rate of 0.62 mm year −1 throughout the catchment. Slight variations are caused by (1) buildup of sediment covers that reduce the weathering rates according to E = E 0 e −S/L where S is sediment thickness and L = 1 m (see Andersen et al., 2015), and (2) the fact that ice cover varies through time in response to the climate forcing. We explored this model to assess the potential of different source areas in explaining the observed trend in 10 Be concentrations (section 4.4.2). We explored three more models to assess the sensitivity of our sampling sites to potential spatial and/or temporal differences in erosion rates (section 4.4.3). Clearly, there are numerous scenarios of how erosion rates may vary in space and time, and it is not the scope of this paper to provide a unique spatiotemporal pattern of erosion rates in the Chhota Shigri catchment. Instead, we refer to established erosion models and explore the range of 10 Be concentrations that they yield during different times of the simulation.
In the second model (Slope), sediment is produced as a function of local slope according to (e.g., Culling, 1960): where K s is a constant. This formulation is usually applied to sediment flux on soil-mantled hillslopes, whereas we are dealing mostly with bedrock hillslopes. However, we are unaware of any accepted law for sediment production on bedrock hillslopes. The purpose here is to simply induce meaningful spatial variability in erosion rates, without claiming the truth. The resulting pattern of erosion rates displays more variability compared to the first model. In the third model (Frost AT), sediment is produced by frost cracking (e.g., Walder & Hallet, 1985), which we computed with the model of Andersen et al. (2015), and using the atmospheric mean annual temperature as the relevant forcing. Consequently, spatial differences in erosion rates are dominated by a vertical gradient (Figure 9c). In the fourth model (Frost LST), sediment is again produced by frost cracking, but this time, the temperature to calculate frost-cracking is derived from mean annual land surface temperatures based on the ASTER AG100 v003 global emissivity and land surface temperature data set (Hulley & Hook, 2011, which has a spatial resolution of 100 m.

Comparing Model and Sample Results
To compare the measured 10 Be concentrations from the samples with the modeled 10 Be concentrations, we had to adjust the locations where the samples were taken within the model. This was necessary because the modeled medial moraine did not stay fixed in space during the model runs, but migrated laterally in response to varying ice discharges from the different tributaries of the Chhota Shigri Glacier. After the modeling, we thus fixed the elevation of the sampling spots on the glacier by means of a contour line and updated their transverse position along that line in order to center them on the modeled medial moraine. We achieved this by quantifying the particle density along the line, and moving the sampling spot to the highest particle concentration multiplied by a penalty function π d ð Þ ¼ e −d=D* , where d is the distance from the relocated sampling position during the last time step, and D * = 100 m is a length scale. Places farther away than 100 m were not considered, as the medial moraine did not migrate by much in between modeling time steps.
For direct comparison of the modeled and measured 10 Be concentrations, we collected all supraglacial particles that are within 100 m of the previously adjusted sample locations, and at selected times during the model run ±100 years. We applied these rather generous bounds to account for uncertainties in sampling locations, times, and particle trajectories. We targeted the modeled erosion rates to be similar to what the

10.1029/2020JF005586
Journal of Geophysical Research: Earth Surface average 10 Be concentration of the samples suggests. However, the different sediment production models induce spatial and temporal variability in erosion rates leading to modeled 10 Be concentrations that are both higher and lower than the measured concentrations. Because we are mostly interested in the relative differences in 10 Be concentrations along the medial moraine, we normalized both the measured and the modeled concentrations at the sampling sites by their average concentration before comparing them.

Numerical Model Results
In the following, we will first present the model results from the four different erosion rate scenarios at 2,000 years model time (Figures 9 and 10). At this time, the modeled glacier reached a steady state ice volume and had an extent approximately equal to the present glacier. These scenarios serve as a backdrop for assessing variability in 10 Be concentrations that may arise from spatial variations in erosion rate due to the different erosion scenarios. We then present model results in which we varied temperatures to assess how the source areas of our samples might change during glacier advance and retreat ( Figure 12). For that purpose, we focus on the uniform erosion rate scenario (model Uniform) to better separate the effects of spatial variations in erosion rate and changing source areas. Finally, we compare our sample results with the combined effects of different erosion scenarios and glacier changes (Figures 13 and 14). Movies of the four simulations are available in Scherler and Egholm (2020).

Spatial Variations in Erosion Rates
The four erosion scenarios introduced in section 5.3 give rise to spatial differences in erosion rate (Figure 9a) that are also reflected in the 10 Be concentration of supraglacial debris (Figure 9b). Compared to the model Uniform, the models Slope and Frost AT yield higher erosion rates on hillslopes surrounding the accumulation zone of the glacier, i.e., in the source areas of supraglacial debris. The spatial distribution of land surface temperatures in the model Frost LST results in very high erosion rates on some north-facing hillslopes in the eastern part of the catchment, but otherwise erosion rates are mostly lower than in the other models. As a result of the spatial variations in erosion rates, 10 Be concentration of particles at the ice and land surfaces differ substantially between the different models ( Figure 9b). Note also that differences in debris supply result in greater amounts of debris in and on the glacier in the models Slope and Frost AT, and consequently a slightly larger glacier, due to the melt-lowering effect of the debris cover (Equation 5).
To understand how these spatial differences in erosion rates affect 10 Be concentrations of the medial moraines of the glacier, it is instructive to examine particle properties within a cross section of the glacier. Figure 10 provides an oblique view of particles within and on the glacier from model Uniform together with a cross section through the glacier near sample CSM1, and Figure 11 shows the relative differences in particle 10 Be concentrations along this cross section in the four different models. Depending on the position of the source area with respect to the glacier, the particles take different trajectories through the ice. For example, particles eroding from headwalls high up above the southeastern edge of the glacier are deeply buried in the ice and transported close to the bed (yellow points in Figure 10a). These particles have long transport distances and source elevations, and are found at~600 m distance along the cross section, at a depth of >200 m within the ice (Figures 10b-10e). In contrast, particles from source areas that are close to the cross section, like the source area we infer for our samples (Figure 4), stay at or close to the surface (blue points in Figure 10a).

Journal of Geophysical Research: Earth Surface
Medial moraines in the model (dark blue particle colors in the ablation area of the glacier in Figure 10a) originate at the confluences of individual glacier tributaries, very similar to the actual glacier (cf., Figure 4). The main medial moraine between the eastern and western branches of the glacier, from which we collected our samples, results from intersection of a thin debris band with the ice surface at the elevation of the cross section (black triangle in Figures 10b-10e). At greater depth within the ice, the debris band widens and occupies a considerable fraction of the glacier's width. Within this debris band, we observe lateral differences in travel distance and transport time, but due to the similar elevations and thus production rates in the source areas, the range of particle 10 Be concentrations in the debris band is rather low (Figure 11, model Uniform). We observe greater diversity in particle source elevations in the adjacent medial moraines (Figures 10b-10e), especially those coming from the eastern branch of the glacier, and this diversity goes along with a greater range in particle 10 Be concentrations (Figure 11, model Uniform). The relative differences in particle 10 Be concentrations within the cross section increase with spatially more variable erosion rates ( Figure 11). In particular, small patches of high or low erosion rates are often juxtaposed in the model Slope and result in large spatial gradients in particle 10 Be concentrations. In contrast, erosion rates vary rather smoothly with elevation in the model Frost AT and thus spatial gradients in 10 Be concentrations within the cross section follow mostly the source elevation (cf., Figure 10e).

Effect of Glacier Response to Climate Change
In the following, we will first present the general behavior of the modeled glacier, and subsequently focus on the medial moraine and the 10 Be concentrations. Figure 12 shows the ice thickness and the distribution of supraglacial particles during these times for model Uniform. Compared to the actual glacier (see Figure 4), the modeled glacier at 2,000 and 3,000 years model time has similar extent and supraglacial debris cover. At 2,500 years, the glacier reaches a larger extent and the relative amount of supraglacial debris cover has

10.1029/2020JF005586
Journal of Geophysical Research: Earth Surface visibly decreased (Figure 12b). After the final warming the glacier attains its minimum extent at 3,500 years, whereas the debris cover appears to attain its maximum relative extent.
Despite large variations in glacier extent during the simulation, the range of 10 Be concentrations at the sampling points along the main medial moraine is limited to~4-8 × 10 4 atoms g −1 , with a median concentration close to~6 × 10 4 atoms g −1 (Figure 13, top row). The likely reason for this narrow and relatively static range of 10 Be concentrations is the limited altitudinal extent of the source area that feeds the main medial moraine (Figures 12b and 3; cf., section 4.4.2). This pattern is different at other places on the glacier surface. For example, medial moraines to the east of the one we sampled are comprised of particles that are derived from source areas that span larger elevation ranges, which results in a wider range of 10 Be concentrations (Figure 9). The biggest change in source areas of our samples occurs toward the end of the simulation, when the glacier has degraded enough to allow our lowermost sample to intersect with the glacier boundary (Figure 12b, 3,500 years). At this stage the sample incorporates particles from hillslopes west of the glacier, which also carry particles that were derived from ice higher up on the same hillslope. Additionally, the main source areas from previous times (Figure 12b, ≤3,000 years) appear to contribute relatively less particles to our samples, whereas additional sources higher up along the southern margin of the glacier are now also contributing.
The above-described observations from the model with uniform erosion rates apply in a similar, but not identical manner to the results from the models with spatially and/or temporally varying erosion rates. Be concentrations are generally small. The 5 to 95 percentile range in 10 Be concentrations is narrow (~4-7 × 10 4 atoms g −1 ), except for the lowest sample (CSM5) toward the end of the model (>3,000 years) as described above. This pattern changes when introducing spatially more variable erosion rates in models 2-4. In the scenario of slope-dependent erosion rates (Slope), the majority of the particles have similar 10 Be concentrations, but some source areas exist that deliver particles with decisively higher concentrations (~9-10 × 10 4 atoms g −1 ). The impact on the mean 10 Be concentration at the sample locations, however, is limited, but increases toward the end of the model and for lower sample locations ( Figure 13). Sediment production under the influence of frost-cracking processes further increases the spread of particle 10 Be concentrations at the sample locations and also introduces more temporal variability.

Observed Versus Modeled Downglacier Changes in 10 Be Concentration
Our modeling results show that it is difficult to explain the observed downglacier increase in 10 Be concentrations through spatial variability in erosion rates or through transient effects of the glacier's response to climate change and associated changes in source areas; given the erosion rate scenarios that we tested ( Figure 14). Furthermore, the temporal changes in erosion rates that result from the imposed temperature changes, in combination with the adopted frost-cracking law, exert only small influences on 10 Be concentrations. For example, the imposed cooling after 2,000 years in the model Frost AT results in a 37% increase in erosion rates in the entire catchment of the Chhota Shigri Glacier, but in a 20% decrease in erosion rates along headwalls above the accumulation area. Hence, the cooling accounts for only a slight increase in 10 Be concentrations, but in a decrease in 10 Be concentrations when it was getting warmer in the model. In the model Frost LST, erosion rates along headwalls above the accumulation area increased by almost 100% during the cooling, but only by 8% in the source areas of the sampled medial moraine.

Journal of Geophysical Research: Earth Surface
Although the model with slope-dependent erosion rates appears to exhibit a tendency for increasing 10 Be concentrations downglacier, the magnitude of this increase is smaller than the variance of 10 Be concentrations among particles at each sample location ( Figure 14). The models in which erosion rates are driven by frost cracking on the other hand-if anything-rather appear to exhibit a tendency for decreasing 10 Be concentrations. But these models, too, display spatial variability in 10 Be concentrations that is smaller than at-a-sample variance. We observed the greatest changes in 10 Be concentrations at the sample locations when particles with distinctly different source areas reached the glacier surface. This effect is more pronounced toward the end of the simulation and for the lowermost sample ( Figure 13). We emphasize, however, that we collected all of our samples, including the lowest one, from areas where the sampled medial moraine was clearly separated from the debris of lateral moraines. We thus think that incorporation of debris from source areas other than the inferred one, in between the two main branches of the glacier (Figure 4), is insignificant. As a result, none of our models, at any time during the simulations, yielded the observed downglacier increase in 10 Be concentrations at the sample locations ( Figure 14).
What other mechanisms could potentially explain the observed downglacier increase in 10 Be concentrations? In our modeling, we made two simplifying assumptions about the 10 Be concentration of bedrock and sediment that may actually turn out to be important. The first one is that the surface 10 Be concentration of eroding bedrock is always in steady state with the local erosion rate. This assumption neglects any time-dependent changes of surface 10 Be concentrations due to changes in both erosion rates and ice cover. Such effects could be particularly important during periods of glacier advance and retreat and the associated cover or exposure of bedrock. The surface 10 Be concentration of recently exposed bedrock is controlled by its exposure history to cosmic rays, which is strongly affected by the duration and thickness of the former ice cover, and the potential loss of material due to erosion. Because shielding a bedrock surface with only 10 m of ice already reduces the 10 Be production rate by~99%, it does not require much ice to impose a significant reduction in production rates. Furthermore, because mountain glaciers like Chhota Shigri are most likely wet-based and erode their bed, any inherited 10 Be concentration would be gradually removed.
Assuming the case that the extent of the Chhota Shigri Glacier has gradually declined throughout the Holocene, it is possible that recently exposed bedrock has a significantly lower 10 Be concentration compared to most hillslopes in the source areas of the debris. If, instead, the Chhota Shigri Glacier had a smaller extent prior to the LIA, the material eroded and transported by the glacier during the LIA advance may have had a significantly higher 10 Be concentration due to the erosion of, e.g., lateral moraines and previously more stable hillslopes, compared to later times. In both cases, gradual retreat and LIA advance, we may anticipate that 10 Be concentrations would decrease with time, but the question is at which rate and magnitude? Figure 15 shows longitudinal and transverse sections through the main medial moraine from model Uniform. At 2,000 and 3,000 years model time, the particles at the sample locations have ages betweeñ 200 and 400 years. In other words, if the modeled transport times are correct, our samples were probably eroded from hillslopes between~1600 CE and~1800 CE, but likely somewhat later, because of faster transport rates, and hence younger ages, when the ice was more extensive (see panel 2,500 years in Figure 15). This timing coincides with glacier advances in western Europe (e.g., Holzhauser et al., 2005), but apparently postdates the main phase of LIA-moraine building in the Himalaya (Rowan, 2017). Although such timing appears to be in conflict with the LIA-advance scenario, the recent history of the Chhota Shigri Glacier is too poorly constrained to refute this hypothesis.
In the gradual retreat scenario, thinning of the ice exposes formerly glaciated bedrock, illustrated in Figure 16 for the glacier retreat between 2,500 and 3,000 years in the model Uniform. Although the largest deglaciated areas lie downstream of the retreated glacier terminus, thinning of the ice also exposes formerly glaciated bedrock in source areas of the main medial moraine. These areas would likely contribute material of low 10 Be concentration to the medial moraine, and increasingly so, with continued lowering of the ice surface. Nevertheless, the total amount of such low-concentrated material is likely small, due to the small extent of these areas. Although, repeat Lidar observations from the Kitzsteinhorn in Austria suggest that formerly glaciated bedrock, adjacent to the Bergschrund, may undergo a period of accelerated erosion upon recent exposure (Hartmeyer et al., 2020), thus accounting for a relatively large supply of material, despite the small area. It should be noted that a relative increase in the supply of material from close to the glacier surface would reduce the aggregated 10 Be concentration of debris derived from that hillslope, even without considering former ice cover shielding, simply because of the lower production rates compared to higher up on the hillslope. Another way that material with low 10 Be concentrations may start to contribute to the medial moraine during periods of glacier retreat arises when tributary glaciers detach from the main ice trunk. The associated disruption in particle paths allows subglacial sediment of the tributary glacier to end up on the surface of the glacier that it formerly joined. At present, however, we found no obvious candidate in the main source areas of our samples, where this mechanism might play a significant role.  The second simplifying assumption is that the material eroded from bedrock has a 10 Be concentration that is equal to the surface concentration. This assumption would be unproblematic if surface erosion occurred at a grain-scale. However, rockfalls, rockslides, and debris avalanches typically scour significantly deeper into the bedrock and thus transport material with very different 10 Be concentrations, including material from greater depth with much lower 10 Be concentrations compared to the surface. Calculations by Ward and Anderson (2011) have shown that mismatches between surface concentrations and those obtained by amalgamating many different rock fall clasts are rather low for a wide range of common rock fall size-frequency distributions (Sarr et al., 2019), but they can nevertheless change when rockfall distributions change. Furthermore, the adjustment of bedrock 10 Be concentrations to changes in erosion rate (e.g., due to climatic changes) requires some time. If erosion occurs steadily at the grain scale, the time required to reach isotopic steady state is equal to the time it takes to erode one e-folding length, which is~60 cm in granitic bedrock, and thus~600 years if erosion rates are~1 mm year −1 (Lal, 1991). This timescale can get much longer when erosion occurs sporadically (Parker & Perg, 2005), as can be expected in glacial landscapes (Ward & Anderson, 2011). We therefore do not expect that the observed downglacier trend in 10 Be concentrations could reliably reflect a temporal change, or shift, in erosion rate. However, this reasoning applies only to sustained changes in erosion rate. It is also conceivable that warming temperatures in glacial landscapes trigger a pulse of greatly enhanced erosion, as observed at the Kitzsteinhorn in Austria (Hartmeyer et al., 2020). In fact, warming-related glacier thinning and the degradation of bedrock permafrost in particular are considered to be important controls on recent rock falls and slope failures in many mountain ranges on Earth (e.g., Gruber et al., 2004;Gruber & Haeberli, 2007;Haeberli et al., 1997).
In summary, although we were not able to reproduce the observed trend in 10 Be concentrations in our numerical simulations, the modeling allowed us to exclude several of the processes mentioned in section 4.4. Specifically, we were able to show that spatial variability of 10 Be production rates is too low in the source areas of the main medial moraine (Figure 4) to explain the observed downglacier trend ( Figure 5). In addition, expansion of source areas to parts of the catchment with significantly different 10 Be production rates or erosion rates is unlikely, even during a period of glacier retreat. The different erosion rate scenarios that we explored (Figure 9) are certainly not exhaustive; yet our modeling has shown that none of them yielded substantial spatial and/or temporal changes in erosion rates. However, a pulse of enhanced erosion due to glacier recession and destabilization of hillslopes, possibly from recently deglaciated areas, may be able to explain the observed downglacier trend in 10 Be concentrations. To test this further, our model would require capabilities to track transiently evolving 10 Be concentrations in and on hillslopes as well as stochastic erosion, which is beyond the scope of this paper.

Dynamic Evolution of Supraglacial Debris Cover
Our model results highlight tight linkages of how changes in climate and the associated impacts on glacier mass balance affect the spatial distribution and composition of debris cover. At the scale of the entire glacier, we observed contrasting patterns of changes in ice and debris volume in the four models ( Figure 17). In general, the steady state configuration of supraglacial debris is controlled by the relative amounts of debris and ice that is transported by the glacier, and by the speed at which the glacier transports this admixture (e.g., Anderson & Anderson, 2016;Scherler et al., 2011). Because the climate forcing was the same in all models, differences between the models can be largely linked to differences in the supply of debris. For example, after the initial build-up of debris, climate cooling and ice advance resulted in an increase also in the volume of debris ( Figure 17) for all models, although more marked in models with variable compared to uniform erosion rates. The reason for this is probably twofold. First, an overall higher erosion rate delivers more sediment to the glacier and makes it more susceptible to melt rate-influencing effects of the debris cover. This most likely accounts for the differences between the models Uniform and Slope. Second, in the models Frost AT and Frost LST, erosion rates and sediment delivery are also influenced by climate change, which may thus account for a more varied response of the debris cover during the simulations.
In all models, changes between different steady state ice volumes are rather rapid and followed by prolonged time periods of changes or fluctuations in supraglacial debris cover. Inspection of the model results suggests that these changes in debris cover appear in part to be related to the build-up of debris near the edge of the glacier. Because the ice velocity near the edge is very low, these particles tend to stay on the ice for a long time, potentially building up thick debris cover that reduces surface ablation rates (cf., Equation 5). In consequence, the ice can get thicker, accelerate, and advance one grid cell or more, carrying the debris along. Subsequent thinning may eventually result in loss of a significant portion of the debris particles. While these processes appear to correctly represent ice dynamics, any additional particle movement that is independent from the ice, e.g., by sliding down a steep surface, is not included in the model. In other words, the spatial distribution of supraglacial debris thickness is entirely related to its passive transport by the ice, and the spatial scale at which debris thickness can possibly vary is confined to the pixel size (25 m). In a previous version of the model, we eliminated very old particles at the glacier edge and the resulting changes in debris volume during the simulation were much smoother. Therefore, the variations in debris volume that we observe at steady state ice volumes should be interpreted with caution.
An interesting development in all four models is that during the final warming period (after 3,000 years), the ice volume shrinks substantially, but the debris volume does not decrease by much. In some models, the debris volume appears to be higher than early during the simulation. Such a development reflects the expansion of supraglacial debris cover during warming periods and is consistent with historical observations (e.g., Deline, 2005;Kirkbride, 1993), and with the negative scaling of fractional debris cover with glacier size observed by Scherler et al. (2018) at a global scale. In the case of the Chhota Shigri Glacier, the melt-reducing effect of thick debris cover for the mass balance of the entire glacier is relatively unimportant when the glacier is large, but it becomes more important when the glacier is small, after~3,500 years model time. It is noteworthy that the median particle age throughout most of the simulations is relatively constant, despite the marked changes in debris volumes. This reflects the combined changes in ice velocity and glacier size that determine the residence time of the particles. After~3,500 years, however, median ages of both englacial and supraglacial debris are somewhat older. We expect that this trend continues when warming continues, because the same amount of debris is transported with less ice, which results in high englacial debris concentrations that lead to more debris cover, which lowers surface melt rates and allows the glacier to remain longer, hence increasing the residence time of the debris (cf., Anderson & Anderson, 2016;Banerjee & Shankar, 2013).
At the scale of an individual medial moraine, the composition and structure of the exposed debris can change considerably. To illustrate this point, consider the particle transport times shown in longitudinal and transverse sections through the main medial moraine from the model Uniform at different times (Figure 15). At any given sample location (black triangles), the mean age of supraglacial particles changes with the different glacial extents (model times). Whereas at 2,000 and 3,000 years, the supraglacial particles of the main medial moraine (Figure 15a) have ages between~150 and 300 years, they are onlỹ 100-200 years old at 2,500 years model time, but~300-600 years old at 3,500 years model time. This phenomenon is true also for the other medial moraines (Figure 15b), and is related to the different transport times through the ice, and to the connectivity of tributaries with the trunk glacier. Specifically, tributaries that are only weakly connected to the trunk glacier have little ice discharge and transport debris with long residence times on the tributary glacier. This leads to debris inputs with high 10 Be concentrations, due to slow transport and high production rates at the typically high and steep headwalls of tributary glaciers (Figures 9b and 15).
As a result of these transient changes, the supraglacial debris that one would sample from the ice surface can stem from distinctly different times in the past. During glacier advance, the time lag between debris generation on headwalls and its supraglacial exposure is short, due to fast transport rates. During phases of glacier retreat, this time increases. Moreover, the associated changes in ablation area exposes debris from different origin, potentially leading to changes in 10 Be concentrations that are not necessarily associated to changes in erosion rate. In the case of the Chhota Shigri Glacier medial moraine, however, the changes in source area are small and thus the temporal evolution in erosion rates, due to the different sampling positions, is probably less than~100-150 years ( Figure 15). During this time, the glacier likely had a negative annual mass balance (Engelhardt et al., 2017). If it is true that the observed downglacier increase in 10 Be concentrations at the Chhota Shigri Glacier is related to thinning and exposing formerly ice-covered hillslope areas, we would expect to observe similar increases at other glaciers, where it can be shown that the source area of the exposed debris did not change by much. The time-transgressive nature of debris exposed in medial moraines thus provides an opportunity to reconstruct erosion rate histories, provided that the source areas and the additional 10 Be concentration during transport can be constrained.

Implications for Cosmogenic Nuclide Dating of Moraines
Surface exposure dating of latero-frontal moraines with cosmogenic nuclides has become an important tool for reconstructing paleoclimate in formerly glaciated mountainous regions (e.g., Ivy-Ochs & Briner, 2014;Owen et al., 2009;Phillips et al., 1996). Our model results provide insight into the exposure history of debris that constitutes such moraines (cf., Ward & Anderson, 2011). Specifically, in all our models, debris is typically exposed at the ice surface for several hundred years before it leaves the ice and potentially gets incorporated into a moraine. A large fraction of this time is contained in the final departure from the ice surface close to its margin. But even farther away from the margin, the age of supraglacial debris can be quite high. For our samples collected from the Chhota Shigri Glacier surface, modeled mean debris ages are typically between~200 and 500 years, but they can get as old as~1,000 years for the lowermost sample at 3,500 years model time. Although we did not collect any samples from the latero-frontal moraines of the Chhota Shigri Glacier, we would expect to obtain ages that are on average a few hundred to a thousand years older than the climate change that might have induced the glacier retreat, provided samples of similar grain size as our medial moraine samples.

Methodological Conclusions
Cosmogenic nuclides are nowadays routinely used to estimate catchment-averaged erosion rates from river sediments. Because glaciers can be considered rivers of ice (Elsa of Arendelle, 2019), supraglacial debris can serve a similar purpose. However, both the sampling and the calculation of erosion rates from measured concentrations are less straight forward and require additional considerations. To avoid issues related to the stochastic nature of rock falls and the limited mixing of debris during glacial transport, Ward and Anderson (2011) used the average 10 Be concentration from their samples from a given glacier to estimate an erosion rate. Such an approach may be preferred, if it can be shown that all the samples share the same exposure history, and differences in 10 Be concentration result purely from stochastic differences in exposure time and/or depth on a hillslope. However, if it is true that the current glacier thinning and retreat significantly affect 10 Be concentrations of supraglacial debris, this approach may introduce bias that does not correctly reflect erosion rates in the long-term. In that case, debris that was eroded prior to the onset of glacier thinning may be less affected compared to debris that was recently eroded. However, depending on how much and how fast erosion rates in glacial landscapes change in response to climatic changes, it may also be that isotopic steady state is hardly every reached. More studies of spatial variations in 10 Be concentrations along medial moraines are needed to assess these issues.
The modeling results from our study underscore the need to constrain the transport history of the debris and its source area as well as possible. This is relatively easy for well-defined medial moraines, especially on small glaciers, but much harder for large glaciers with continuous debris blankets. In that case, issues can arise due to unknown source areas that exhibit a wide range of production rates, lateral mixing of supraglacial debris by differential melting due to heterogeneous debris thickness, or complex exposure histories that include periods of temporary burial in supraglacial ponds or hollows. Additional issues can arise if different glacier branches merge in the ablation area and contribute debris to a medial moraine that may have been exposed in lateral moraines for unknown time periods. Moreover, if the potential source areas include tributary glaciers that have detached from the main valley glacier, there exists the risk that material ends up on the glacier surface that has been eroded from subglacial sources. Unfortunately, such complex exposure histories are difficult to constrain quantitatively, and thus to account for when interpreting such samples. At the moment, our model and the new particle tracing module is only partly able to simulate more complex exposure histories of particles. Although we did not consider subglacial erosion and sediment transport in our models, it can be included to account for more complex glacial systems. The mixing of particles within supraglacial debris cover is more difficult and we currently have no representation of surface processes that could move particles around on the ice or crevasses where particles could fall into. However, for a glacier like Chhota Shigri, surface transport and crevasses most likely play a minor role.
Finally, our study has shown the potential that the current shrinking of glaciers imprints itself in cosmogenic nuclide concentrations of supraglacial debris. Although this could render the determination of reliable headwall erosion rates more difficult, it may also allow us to exploit medial moraines as archives and better constrain the response of glacial landscapes to climate change. To make progress along these lines, two approaches appear promising. First, measuring in situ-produced 14 C in medial moraine samples could provide clues on transient adjustments on a shorter time scale than 10 Be. Because of the much shorter half-life of 14 C, the ratio of 14 C to 10 Be is a sensitive indicator of both complex exposure histories due to changes in ice cover and rapid surface erosion events (Hippe, 2017). Second, incorporating the effects of transient changes in erosion rates and ice cover as well as stochastic erosion on cosmogenic nuclide concentrations in the numerical model would allow to better explore the impact of past and present climate change, and to predict, e.g., the statistical variation in nuclide concentration of the debris cover given certain climatic boundary conditions. However, it should be noted that such model capabilities come along with the requirement to define a history for poorly constrained processes such as subglacial erosion or rock fall erosion in order to correctly initialize bedrock nuclide concentrations. In summary, there exist ample opportunities to study the response of glacial landscapes to climate change by coupling cosmogenic nuclides with numerical models, and we hope that our work promotes further studies on other glaciers.

10
Be concentrations of five amalgamated debris samples from the main medial moraine of the Chhota Shigri Glacier, Indian Himalaya, suggest headwall erosion rates that are~0.5-1 mm year −1 . Apart from the lowermost sample, the 10 Be concentrations increase downglacier by values that cannot be explained through simple passive transport on the ice. We developed a modeling framework, in which Lagrangian particle tracing allowed us to study the evolution of 10 Be concentrations in glacial debris during transport. Although we were not able to reproduce the observed downglacier increase in 10 Be concentrations in our numerical experiments, the modeling yielded valuable insights into the processes that control spatial and temporal variability in 10 Be concentrations of supraglacial debris. In particular, the modeling has shown that large changes in 10 Be concentration can be expected in places that expose debris from source areas that cover wide ranges in elevation and/or erosion rate. However, this can be excluded for the main medial moraine of the Chhota Shigri Glacier. At the expense of alternative explanations, we are currently inclined to suggest that the observed downglacier increase in 10 Be concentrations is due to exposure and erosion of recently deglaciated parts of hillslopes. To account for this explanation, future coupled modeling of ice, debris, and cosmogenic nuclides should thus include transient effects on bedrock 10 Be concentrations. All of the discussed reasons for potential downglacier changes in the 10 Be concentration are not limited to our samples or the Chhota Shigri Glacier. In fact, our modeling shows that the 10 Be concentration of debris exposed in medial moraines of any valley glacier is a complex function of its source area, its transport path, and any transient responses of hillslopes and glaciers to perturbations of the steady state.

Data Availability Statement
Landsat 8 imagery provided by courtesy of the U.S. Geological Survey. The source code and input files of the numerical model, results from the numerical experiments, and the sample data are available from Scherler and Egholm (2020) and Egholm and Scherler (2020).