# Modelling the effects of ice transport and sediment sources on the form of detrital thermochronological age probability distributions from glacial settings

Maxime Bernard, Philippe Steer, Kerry Gallagher, David Lundbek Egholm

## ► To cite this version:

Maxime Bernard, Philippe Steer, Kerry Gallagher, David Lundbek Egholm. Modelling the effects of ice transport and sediment sources on the form of detrital thermochronological age probability distributions from glacial settings. Earth Surface Dynamics, 2020, 8 (4), pp.931-953. 10.5194/esurf-8-931-2020. insu-02998746

HAL Id: insu-02998746

https://insu.hal.science/insu-02998746

Submitted on 10 Nov 2020

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Modelling the effects of ice transport and sediment sources on the form of detrital thermochronological age probability distributions from glacial settings

Maxime Bernard¹, Philippe Steer¹, Kerry Gallagher¹, and David Lundbek Egholm²

¹Univ Rennes, CNRS, Géosciences Rennes, UMR 6118, 35000 Rennes, France

²Department of Geoscience, Aarhus University, Aarhus, Denmark

**Correspondence:** Maxime Bernard (mabernard.ns@gmail.com)

Received: 4 February 2020 – Discussion started: 11 February 2020

Revised: 4 June 2020 – Accepted: 6 September 2020 – Published: 10 November 2020

**Abstract.** The impact of glaciers on the Quaternary evolution of mountainous landscapes remains controversial. Although in situ or bedrock low-temperature thermochronology offers insights on past rock exhumation and landscape erosion, the method also suffers from potential biases due to the difficulty of sampling bedrock buried under glaciers. Detrital thermochronology attempts to overcome this issue by sampling sediments at e.g. the catchment outlet, a component of which may originate from beneath the ice. However, detrital age distributions not only reflect the catchment exhumation, but also spatially variable patterns and rates of surface erosion and sediment transport. In this study, we use a new version of a glacial landscape evolution model, iSOSIA, to address the effect of erosion and sediment transport by ice on the form of synthetic detrital age distributions. Sediments are tracked as Lagrangian particles formed by bedrock erosion, and their transport is restricted to ice or hillslope processes, neglecting subglacial hydrology, until they are deposited. We base our model on the Tiedemann Glacier (British Columbia, Canada), which has simple morphological characteristics, such as a linear form and no connectivity to large tributary glaciers. Synthetic detrital age distributions are generated by specifying an erosion history, then sampling sediment particles at the frontal moraine of the modelled glacier. Results show that sediment sources, reflecting different processes such as glacier and hillslope erosion, can have distinct bedrock age distribution signatures, and estimating such distributions should help to identify predominant sources in the sampling site. However, discrepancies between the detrital and bedrock age distributions occur due to (i) the selective storage of a large proportion of sediments in small tributary glaciers and in lateral moraines, (ii) the large range of particle transport times due to varying transport lengths and strong variability of glacier ice velocity, (iii) the heterogeneous pattern of erosion, and (iv) the advective nature of glacier sediment transport along ice streamlines. This last factor leads to a poor lateral mixing of particle detrital signatures inside the frontal moraine, and then local sampling of the frontal moraine is likely to reflect local sources upstream. Therefore, sampling randomly across the moraine is preferred for a more representative view of the catchment age distribution. Finally, systematic comparisons between synthetic (U–Th)/He and fission track detrital ages, with different bedrock age-elevation profiles and different relative age uncertainties, show that the nature of the age-elevation relationship and age uncertainties largely control the ability to track sediment sources in the detrital record. However, depending on the erosion pattern spatially, qualitative first-order information may still be extracted from a thermochronological system with high uncertainties (> 30 %). Overall, our results demonstrate that detrital age distributions in glaciated catchments are strongly impacted not only by erosion and exhumation but also by sediment transport processes and their spatial variability. However, when combined with bedrock age distributions, detrital thermochronology offers a novel means to constrain the transport pattern and time of sediment particles.
# 1 Introduction

Glaciers have left a profound impact on the topography of mountainous landscapes, in particular by eroding deep glacial valleys and depositing a large volume of sediments in moraines. As glacier dynamics are linked to climate, an active area of research aims to characterize the role of glacial erosion on the dynamics and relief development of mountain belts during the recent Quaternary glaciations (e.g. Zachos et al., 2001; Molnar and England, 1990; Beaumont et al., 1992; Montgomery and Brandon, 2002; Brozović et al., 1997; Whipple, 2009; Steer et al., 2012; Champagnac et al., 2014). To address these questions, two timescales have typically been considered: a longer timescale ($10^5–10^6$ years) to assess the potential glacial imprint on the landscape and a shorter timescale ($10^1–10^4$ years) to understand how ice actually erodes the landscape. For example, some studies integrated glacial sediment records worldwide and in situ low-temperature thermochronology data to estimate glacial erosion rates. They showed average erosion rates of $10^0–10^3$ mm yr$^{-1}$ on short timescales ($10^3–10^5$ years) and long-term ($> 10^6$ years) average erosion rates of $10^{-2}–10^0$ mm yr$^{-1}$ (Hallet et al., 1996; Koppes and Montgomery, 2009; Valla et al., 2011b; Koppes et al., 2015; Bernard et al., 2016).

Herman et al. (2013) used a global compilation of in situ thermochronological data and an inverse approach to infer an increase in erosion rates for all mountain ranges in the Quaternary period. They suggested that this effect is more pronounced for glaciated mountains, implying a significant role of glaciation on erosion rates. However, the results of this study have been contested (Willenbring and Jerolmack, 2016; Schildgen et al., 2018) and thus the conclusion debated. More recently, Yanites and Ehlers (2016) correlated high glacier ice-sliding velocities with high denudation rates deduced by thermochronological data in the southern Coast Mountains, British Columbia. They also found that glacial erosion may only occur for less than 20% of a glacial-interglacial cycle in some areas and may explain the discrepancy between the longer and shorter timescale for erosion rates. In situ thermochronology consists of collecting bedrock samples from discrete locations to map thermochronological ages and exhumation rates (e.g. Fitzgerald and Stump, 1997; Valla et al., 2011b; Herman et al., 2013). However, this approach offers a potentially biased assessment of catchment-wide exhumation patterns as it is dependent on the spatial bedrock sampling strategy (Ehlers, 2005; Valla et al., 2011a). For example, the age-elevation relationship inferred from in situ or bedrock thermochronology data may not capture younger ages expected along the valley floor and buried under the ice (Enkelmann et al., 2009; Grabowski et al., 2013).

In glacial environments, erosion mostly occurs subglacialy through abrasion and quarrying (Boulton, 1982; Hallet et al., 1996) or supraglacialy (i.e. ice-free areas) through periglacial mechanisms and gravitational processes (Matsuoka and Murton, 2008). These two sources (i.e. supraglacial and subglacial areas), and their contribution to erosion, define the relative proportion of sediments that enter the glacier transport system (Small, 1987). Sediments are transported by glaciers by (i) subglacial water flow through cavity or channel systems (Kirkbride, 2002; Alley et al., 1997; Spedding, 2000) and (ii) by ice internal deformation for sediments incorporated within or above the ice (Hambrey et al., 1999; Goodsell et al., 2005).

Therefore, sampling and analysing these sediments with detrital thermochronology have the potential to provide a more spatially integrated view of the catchment erosion pattern on short timescales ($10^1–10^4$ years) and thereby a more representative indication of how ice is eroding the landscape (e.g. Stock et al., 2006; Tranel et al., 2011; Ehlers et al., 2015). The detrital sampling protocol focusses on collecting glacial deposits, generally close to the outlet of the catchment (Ruhl and Hodges, 2005; Stock et al., 2006; Falkowski et al., 2016; Glotzbach et al., 2018). A major motivation for detrital sampling is that we can potentially obtain grains from regions inaccessible for bedrock sampling due to ice cover and lack of outcrop or, more pragmatically, logistical considerations (e.g. cost). However, a priori knowledge about the distribution of thermochronological ages with elevation (e.g. Brewer et al., 2003), and the mineral fertility of the sources (Moecher and Samson, 2006), is often required to reliably interpret the data from detrital thermochronology.

Thermochronological age distributions of detrital samples are generally presented as synoptic probability density functions (SPDF) (Brewer et al., 2003; Ruhl and Hodges, 2005). Analysis of many glaciated catchments has been based on the interpretation of these SPDFs (Stock et al., 2006; Tranel et al., 2011; Avdeev et al., 2011; Thomson et al., 2013). More recently, Enkelmann and Ehlers (2015) investigated the degree of mixing of sediments from ice cores at the terminus of the Tiedemann Glacier (Coast Mountains, Canada) using both detrital and in situ low-temperature thermochronology. Comparison of the detrital age distributions from the ice-cored terminal moraine and from glacial outwash (Ehlers et al., 2015) showed that sampling through the terminal moraine is similar to sampling the proglacial river, supporting the idea of an efficient vertical mixing of glaciated sediments at the glacier front. The latter result is also consistent with results obtained on the Malaspina Glacier, Alaska (Enkelmann et al., 2009; Grabowski et al., 2013). However, Enkelmann and Ehlers (2015) pointed out the sensitivity of their results to the relatively high uncertainties on the ages. Furthermore, while the shape of detrital SPDFs is expected to be mainly controlled by the catchment exhumation history,
other processes such as grain erosion, transport, and deposition are likely to play a role. Consequently, characterizing the effect of surface processes on the shape of detrital SPDFs may help improve their interpretation.

In this study, we present a numerical approach that allows us to explore the effect of sediment transport by ice and the role of different source areas (i.e. subglacial and supraglacial) in shaping detrital SPDFs at a glacier front. To this end, we combined a new version of a glacial landscape evolution model, iSOSIA (Egholm et al., 2011), that allows for the tracking of sediments with an external routine (Egholm and Bernard, 2020) that calculates the detrital thermochronological age distributions. We chose to apply our numerical approach to the Tiedemann Glacier catchment, as it shows simple morphological characteristics and because thermochronological data are available (Enkelmann and Ehlers, 2015; Ehlers et al., 2015). We produced our bedrock synthetic thermochronological data from two low-temperature thermochronological systems, apatite (U-Th)/He and apatite fission track (i.e. AHe and AFT, respectively), according to the age-elevation profiles from Enkelmann and Ehlers (2015) for the AFT ages and from Ehlers et al. (2015) for the AHe ages. First, we investigate the kinematics of our sediment transport model in the modelled catchment to assess the role of the ice transport. Next, we focus on the role of the sediment sources by calculating detrital SPDFs resulting from uniform catchment erosion across the catchment as well as hillslope-only sources and glacial-only eroded sources. Finally, we consider non-uniform catchment erosion by letting iSOSIA and its implicit erosion laws control the pattern of erosion. This last experiment allows us to assess the record of complex erosion processes potentially contained in detrital SPDFs. For all experiments, the sediments are transported by ice and hillslopes so that transport by the hydrological system is neglected. The two thermochronological methods (i.e. AFT and AHe) considered have different thermal sensitivities and thus different spatial distributions of bedrock ages in the landscape. They typically have different age uncertainties, allowing us to address the effect of these uncertainties on the interpretation of SPDFs.

# 2 Methods

## 2.1 The glacial landscape evolution model: iSOSIA

The iSOSIA model is a one-layer depth-integrated second-order shallow ice approximation model (Egholm et al., 2011). It is able to simulate ice flowing on relatively steep topography through the computation of membrane stresses (Hindmarsh, 2006) and with a depth-integrated ice flow velocity. These two characteristics make iSOSIA more accurate than models based on the zeroth-order shallow ice approximation, while being computationally efficient compared to models that resolve the full 3D set of Navier-Stokes equa-

tions (Elmer/Ice; Braedstrup et al., 2016). A full description of iSOSIA and the relevant ice flow equations can be found in Egholm et al. (2011, 2012a and b) and Braedstrup et al. (2016).

Glacial hydrology plays an important role in ice-sliding velocity (Clarke, 1987) and has recently been investigated in numerical models (Schoof, 2005, 2010; Iverson, 2012; Ugelvig et al., 2016, 2018). Here we adopt a simplified version of the sliding law of Schoof (2005) that accounts for the opening of cavities due to the roughness of the bed. This defines the area on which the basal shear stress ($\tau_b$) is applied and the basal sliding speed ($u_b$) as

$$u_b = C_s \cdot \frac{\tau_b^n}{C}, \quad (1)$$

where $C_s$ is the ice-sliding constant and $C = 1 - \frac{s}{L}$ defines the proportion of area on which $\tau_b$ is applied, where $L$ is the length between two topographic steps (Fig. S1). The steady-state length of cavities ($S$) is controlled by the effective mass balance (i.e. accounting for basal and internal melting of the ice) that determines the annual average water flux over the glacier, $q_w$. This is specified as follows:

$$S = \left( \frac{L_s \cdot q_w}{k_0 \cdot \sqrt{\nabla \psi}} \right)^{0.8} \cdot \frac{1}{\beta \cdot h_s}, \quad (2)$$

where $L_s$ is the mean cavity spacing in a cell (see Table 1 for all parameter values), $k_0$ is the minimum hydraulic conductivity, $\nabla \psi = \rho_w g \nabla h_{ice}$ is the hydrological head, $\beta$ is a scaling factor, and $h_s$ is the mean height of the topographic step and depends linearly on the slope in the direction of sliding (see Ugelvig et al., 2016; Eq. 10). The term in parentheses on the right-hand side of Eq. (2) represents the cross-sectional area of a cavity ($A_s$). The steady cavity size controls the effective pressure ($N$) of the system and influences the basal sliding speed. However, due to how cavities dictate the volume of water stored at the bed, the effective pressure is also influenced by the opening and closure rate of cavities following

$$N = B \left(\frac{8}{\pi}\right)^{\frac{1}{n}} \left[ \frac{u_b h_s - \beta h_s \frac{\partial S}{\partial t}}{S^2} \right]^{\frac{1}{n}}, \quad (3)$$

where $B$ and $n$ are the ice creep parameters. This hydrological model is similar to that used in Ugelvig et al. (2016), and we refer the reader to this study for more details. We calibrated our hydrological model through a trial-error process by varying $k_0$ (Eq. 2) to produce a reasonable value of basal ice-sliding velocities (i.e. $\times 10 \text{ m yr}^{-1}$; Table S1).

In iSOSIA, the mass balance $M$ is linearly proportional to the mean annual air temperature $T_{air}$:

$$M = \begin{cases} -m_{acc} T_{air}, & \text{if } T_{air} \le 0 \\ -m_{abl} T_{air}, & \text{if } T_{air} > 0' \end{cases} \quad (4)$$
<table><caption>Table 1. List of parameters described in the text.</caption><thead><tr><th>Parameters</th><th>Description</th><th>Value and unity</th></tr></thead><tbody><tr><td colspan="3"><strong>Climate</strong></td></tr><tr><td><i>T</i><sub>sl</sub></td><td>Sea level temperature</td><td>9 °C</td></tr><tr><td><i>d</i><i>T</i><sub>air</sub></td><td>Lapse rate</td><td>0.006 °C m<sup>-1</sup></td></tr><tr><td><i>m</i><sub>acc</sub></td><td>Accumulation gradient</td><td>0.08 m yr<sup>-1</sup> °C<sup>-1</sup></td></tr><tr><td><i>m</i><sub>abl</sub></td><td>Ablation gradient</td><td>0.25 m yr<sup>-1</sup> °C<sup>-1</sup></td></tr><tr><td colspan="3"><strong>Ice</strong></td></tr><tr><td><i>n</i></td><td>Ice creep constant</td><td>3.0</td></tr><tr><td><i>B</i></td><td>Ice viscosity constant</td><td>73.3 × 10<sup>6</sup> Pa s<sup>1/3</sup></td></tr><tr><td><i>C</i><sub>s</sub></td><td>Ice-sliding constant</td><td>1.29 × 10<sup>-2</sup> Pa s<sup>1/3</sup></td></tr><tr><td><i>h</i><sub>ice</sub></td><td>Ice surface elevation</td><td>Model outcome</td></tr><tr><td><i>H</i></td><td>Ice thickness</td><td>Model outcome</td></tr><tr><td>∇<i>b</i><sub>s</sub></td><td>Bed slope in the direction of ice flow</td><td>Dimensionless</td></tr><tr><td><i>k</i><sub>q</sub></td><td>Quarrying coefficient for erosion</td><td>1.97 × 10<sup>-4</sup> s Pa<sup>-3</sup></td></tr><tr><td colspan="3"><strong>Hillslope</strong></td></tr><tr><td><i>S</i><sub>c</sub></td><td>Critical slope</td><td>1.4</td></tr><tr><td><i>K</i><sub>h</sub></td><td>Hillslope diffusivity</td><td>5 m<sup>2</sup> yr<sup>-1</sup></td></tr><tr><td><i>k</i><sub>eh</sub></td><td>Hillslope erosion coefficient</td><td>3.95 × 10<sup>-5</sup> m<sup>2</sup> yr<sup>-1</sup></td></tr><tr><td colspan="3"><strong>Particles</strong></td></tr><tr><td><i>H</i><sub>s</sub></td><td>Particle sediment thickness</td><td>0.01 m</td></tr><tr><td><i>z</i><sub>p</sub></td><td>Burial depth of particles relative to the ice surface</td><td>Dimensionless</td></tr><tr><td colspan="3"><strong>Hydrology</strong></td></tr><tr><td><i>N</i></td><td>Effective pressure</td><td>Model outcome</td></tr><tr><td><i>h</i><sub>s</sub></td><td>Bed step height</td><td>Model outcome</td></tr><tr><td><i>S</i></td><td>Cavity size</td><td>Model outcome</td></tr><tr><td><i>L</i></td><td>Bed step length</td><td>Model outcome</td></tr><tr><td><i>L</i><sub>s</sub></td><td>Cavity spacing</td><td>4 m</td></tr><tr><td><i>k</i><sub>0</sub></td><td>Minimum hydraulic conductivity</td><td>3.10 × 10<sup>-4</sup> kg<sup>-1/2</sup> m<sup>3/2</sup></td></tr><tr><td><i>q</i><sub>w</sub></td><td>Water flux</td><td>Variable</td></tr><tr><td>∇<i>ψ</i></td><td>Hydrological gradient</td><td>∇<i>ψ</i> = <i>q</i><sub>w</sub> ∇<i>h</i><sub>ice</sub></td></tr><tr><td><i>β</i></td><td>Cavity shape parameter</td><td>0.7</td></tr></tbody></table>

where $m_{acc}$ and $m_{abl}$ are the mass accumulation and the ablation gradients with respect to temperature, respectively. $T_{air}$ is assumed to decrease linearly with elevation $h$ with a lapse rate $dT_{air}$ as

$$T_{air} = T_{sl} + dT_{air} \cdot h, \qquad (5)$$

where $T_{sl}$ is the sea level temperature set constant at 9°C. As mentioned in the Introduction, the two main mechanisms of glacial erosion are abrasion and quarrying (Boulton, 1982; Hallet et al., 1996). As debris resulting from abrasion is mostly small in size (i.e. < 50 µm; Hallet, 1979), we follow MacGregor et al. (2009) and consider this debris to be instantaneously transported away from the glacial catchment by meltwater. However, quarrying acting on the lee side of bedrock steps by plucking can produce larger debris. Therefore, in this study we consider sediments to only be pro-

duced by quarrying. Quarrying is generally accounted for by a classical basal-sliding power law, modified to account for some factors directly related to plucking (MacGregor et al., 2000, 2009; Kessler et al., 2008; Egholm et al., 2009; Herman et al., 2011). These factors include pre-existing fractures in the bedrock, effective pressure, the regional bed slope in the direction of ice flow, variation in meltwater drainage, and the ice-bed contact area determined by the opening of cavities (Iverson, 2012; Cohen et al., 2006; Ugelvig et al., 2016, 2018; Anderson, 2014). We follow Ugelvig et al. (2016) in modelling the quarrying rate $E_q$ as

$$E_q = k_q N^3 u_b (\nabla b_s + b_{s0})^2, \qquad (6)$$

where $k_q$ is a scaling coefficient for the effective pressure, $N$, that represents the effect of bedrock lithology and fractures, $\nabla b_s$ is the bed slope in the direction of sliding, and $b_{s0}$ is a
term to allow for negative bed slopes. As we focus on glacier-induced erosion, we ignore fluvial erosion. However, we do consider hillslope erosion, which represents the supraglacial source of sediments. We assume that hillslope erosion, $E_h$, depends non-linearly on the bed gradient, as defined in Andrew and Bucknam (1987) and Roering et al. (1999):

$$E_h = - \frac{k_{eh} \cdot \nabla b}{1 - \left( \frac{|\nabla b_{ij}|}{S_c} \right)^2} \frac{dt}{dl}, \quad (7)$$

where $k_{eh}$ is an erodibility constant, $S_c$ the critical slope, $dt$ the time step, and $dl$ the resolution of a cell ($dl = 100$ m).

## 2.2 Particle tracking and sampling

To calculate the detrital age distributions produced from a glaciated catchment, we need to track sediments from the source to their final site of deposition. Accordingly, we define “particles” as sediment trackers in our models. Sediments are produced as a result of local erosion. When sediment thickness exceeds a threshold $H_s = 0.01$ m, a particle is generated and can be transported away from its initial cell. This threshold is chosen to lead to a reasonable number of particles within models while keeping a sufficient time resolution in the generation of particles (i.e. with an erosion rate of $1 \text{ mm yr}^{-1}$ and $H_s = 0.01$ m, a particle is generated every 10 years in a cell). Depending on its position within the catchment, a particle is transported by different processes. Similar to Eq. (7), we assume that the flux, $q_{ph}$, of an unglaciated particle depends non-linearly on the bed gradient ($\nabla b$):

$$q_{ph} = - \frac{K_h \cdot \nabla b}{1 - \left( \frac{|\nabla b|}{S_c} \right)^2}, \quad (8)$$

where $K_h$ is the hillslope diffusion coefficient. Glaciated particles can be incorporated into the ice either from the ice surface, for supraglacial debris, or from the ice-bed interface by basal quarrying. The horizontal velocity of each glaciated particle is computed as a function of its depth within the ice, $z_p$, relative to the ice surface. Because iSOSIA is a one-layer depth-integrated model, the horizontal velocity depth profile is computed following Glen's flow law (Glen, 1952). We follow Rowan et al. (2015) in assuming that, due to the ice viscosity, the horizontal ice velocity decays as a fourth-order polynomial of the ice thickness:

$$u_p(z_p) = \frac{5}{4} [1 - z_p^4] \bar{u} + u_b, \quad (9)$$

where $\bar{u} = \bar{u}_d + u_b$ is the average velocity of the ice, with $u_b$ the basal ice speed and $\bar{u}_d$ the average velocity due to the internal deformation of the ice, which is approximated as a fourth-order polynomial of the local ice thickness as well as the bed and ice surface gradients (Braedstrup et al., 2014).

We stress again that we neglect the transport of particles by meltwater and focus only on the ice and hillslope processes. Velocity profiles for glaciated and unglaciated particles are depicted in Fig. S2.

The burial depth of particles ($z_p$) within the ice is expressed relative to the ice surface (i.e. $z_p = 0$ at the ice surface and $z_p = 1$ at its base; Fig. 1) and computed according to the mass balance:

$$\frac{\partial z_p}{\partial t} = \frac{z_p \dot{M}_b - (1 - z_p) \dot{M}_s}{H}, \quad (10)$$

where $\frac{\partial z_p}{\partial t}$ is the change in burial depth of a particle through time, $\dot{M}_b$ and $\dot{M}_s$ are the surface and basal melting at time $t$, respectively, and $H$ is the ice thickness. Relative to the ice surface elevation, particles are moved downward according to the accumulation rate (in the accumulation zone $\dot{M}_s < 0$) or upward with the lowering of the ice surface due to melting (in the ablation zone $\dot{M}_s > 0$; Fig. 1).

In the experiments presented in the following sections, detrital particles are sampled throughout the frontal moraine, defined as the glacier tongue (Figs. 1 and 5d). As vertical mixing of sediment has been reported in some glacier fronts (Enkelmann et al., 2009; Grabowski et al., 2013; Enkelmann and Ehlers, 2015), we sample particles independently of their vertical position within the ice thickness to infer detrital age distributions.

## 2.3 Modelling the Tiedemann Glacier

Our main goal is to investigate the influence of sediment sources and transport by ice on the shape of detrital SPDFs at the glacier front. We present our modelling approach using the Tiedemann Glacier catchment (Coast Mountains, British Columbia, Canada; Fig. 2) as a reference case. This region has been previously studied using detrital and in situ thermochronology (Ehlers et al., 2015; Enkelmann and Ehlers, 2015). Moreover, the Tiedemann Glacier has simple morphological characteristics with a linear main glacial valley and small-extent tributary glaciers. We calibrated our iSOSIA model to mimic the Tiedemann Glacier through a trial-error process. We varied climatic, ice, and hydrological parameters ($dT_a$, $C_s$, $M_{acc}$, and $k_0$; see Table S1) and compared the resulting mean ice thickness and location of the glacier tongue against the results of the ITMIX experiments (Farinotti et al., 2017, 2019), which predicts the thickness of many glaciers around the world and thus the Tiedemann Glacier (Fig. S3). To produce synthetic detrital age distributions resulting from this model, we restrict the glacier to steady state and impose a constant topography (i.e. erosion produces particles but no topographic changes). We stress that our goal is not to fit the dynamics of the actual Tiedemann Glacier, which is in a phase of retreat (Tennant et al., 2012). Rather, we ask how surface processes affect the shape of detrital SPDFs in a simple case of constant and steady transport and erosion within a glaciated catchment.
<figure><img src="image_3.png" /></figure>

**Figure 1.** Particle transport by ice. Particles originating from hillslopes enter the glacier system from the surface (z<sub>p</sub> = 0), whereas particles produced by ice erosion are incorporated into the basal ice (z<sub>p</sub> = 1). The transport of glaciated particles follows the ice flow pattern from the sources to the frontal moraine. ELA is the equilibrium line altitude.

<figure><img src="image_4.png" /></figure>

**Figure 2.** Satellite view of the Tiedemann Glacier, British Columbia, Canada. The red contour represents the catchment area. The inset shows the regional location in the Coast Mountains, Canada. Note the dominantly linear pattern of sediment transport (image from © Google Earth, Maxar Technologies; inset map with data from Jarvis et al., 2008).

The elevation of the Tiedemann Glacier catchment ranges from 363 to 3920 m for a total local relief of 3557 m. The resulting maximum ice thickness is 454 m, with a mean ice velocity of 13.2 m yr<sup>-1</sup> and a range of 0 to 109.2 m yr<sup>-1</sup> (Fig. 3), which are characteristic of Alpine glaciers (Benn and Evans, 2014). The mean ice velocity is mainly controlled by the basal sliding speed, which in turn is influenced by the shear stress on the bed and the distribution of the effective pressure (Fig. S4). This results in two main regions of relatively high velocities (~60 m yr<sup>-1</sup>) located at ~7 and ~14 km in the glacier latitudinal coordinates (Fig. 3c).

## 2.4 Synthetic ages and production of detrital age distributions

To produce detrital SPDFs we need to set the bedrock or in situ age distributions. To this end, we consider the age-elevation profiles from two low-temperature systems presented in Enkelmann and Ehlers (2015) (AFT ages) and from Ehlers et al. (2015) (AHe ages; Fig. 4a). We assume that true AHe and AFT ages in the bedrock at a given elevation are given by the age-elevation relationship. In practice, thermochronological analysis on a single grain is presented as an age and its associated uncertainty. Thus, the measured age can be considered a noisy sample of the true age. To produce an SPDF, we simulate this process for the detrital grain ages by adding random noise to the true ages and assume the detrital grain ages are not modified during transport and deposition.
<figure><img src="image_5.png" /><figcaption>Figure 3. Characteristics of the Tiedemann Glacier model. (a) The underlying topography, (b) ice thickness, and (c) mean ice velocity. The mean ice thickness has been calibrated with the ITMIX experiments (see Fig. S3).</figcaption></figure>

In Ehlers et al. (2015), AHe bedrock ages and thus detrital data have relatively high uncertainties (~21%). However, we stress here that one of our goals is to assess the effect of age uncertainty on the interpretation of SPDFs, so we chose to apply lower uncertainty for the synthetic AHe ages. For this reason, we adopt uncertainties of 10% for the synthetic AHe ages sampled from the age-elevation profile. This is similar to the typical reproducibility in aliquots of single-grain ages (e.g. ~6%; Farley and Stocki, 2002).

The AFT data from Enkelmann and Ehlers (2015) show > 30% relative uncertainties. We generate synthetic AFT ages following an approach similar to Gallagher and Parra (2020). We use the external detector method (EDM) age equation (e.g. Hurford and Green, 1983) to calculate the ratio of the spontaneous and induced track density ($\rho_s/\rho_i$) corresponding to a given AFT age in the age-elevation profile. To allow for the effect of variable high uncertainties, we simulate the range of relative error of the Enkelmann and Ehlers (2015) study by randomly sampling the published values of $N_s+N_i$ (i.e. the sum of the number of spontaneous and induced tracks). We then use a binomial distribution with the parameter $\theta = \rho_s/\rho_i$ to sample $N_s$ values, corresponding to the number of spontaneous tracks, given the sampled value of $N_s+N_i$ (see Gallagher, 1995; Galbraith, 2005). With the randomly sampled value of $N_s$, conditional on the sampled value of $N_s+N_i$, we easily obtain a value for $N_i$. We can estimate a “noisy” AFT age and relative error from the sampled $N_s$ and $N_i$ values by using the EDM equation. For both the AHe and

AFT systems, we attribute to each particle a single age and uncertainty in accordance with its source location elevation.

The probability age distributions are then produced by using Gaussian kernel density estimate (Brandon, 1996; Ruhl and Hodges, 2005) with the mean and standard deviation equal to the age and error for both the AHe and AFT ages following

$$SPDF_k = \frac{1}{N_p \alpha \sigma_{t_k} \sqrt{2\pi}} \exp \left[ -\frac{1}{2} \left( \frac{t - t_k}{\alpha \sigma_{t_k}} \right)^2 \right], \quad (11)$$

where $N_p$ is the total number of sampled particles (i.e. ages), $t$ is the range of ages in which we compute the SPDF, $t_k$ is the age of a particle considered, $\sigma_t$ is the associated age uncertainty, and $\alpha = 0.6$ is a scaling factor that handles the resolution and the precision of age components in an SPDF (Brandon, 1996). To compare age distributions, we also define a cumulative synoptic probability density function, which is simply the cumulated sum of the SPDF:

$$CSPDF = \sum_{j=0}^{N_p} SPDF_j. \quad (12)$$

Given the relatively high number of particles within the frontal moraine (i.e. > $10^6$) and to mimic real detrital sampling, all the detrital SPDFs presented in the following sections are produced by randomly sampling 105 particles from the total, independently of their source origin (Sect. 3.2, 3.3,
<figure><img src="image_6.png" /></figure>

**Figure 4.** Synthetic bedrock thermochronological age distributions for the Tiedemann Glacier model. The age-elevation profile for the two thermochronological systems AHe and AFT (**a**). Black and white circles represent the bedrock AHe and AFT ages from Ehlers et al. (2015) and Enkelmann and Ehlers (2015), respectively. Solid black lines are the best-fit solution for the data. Shaded areas for the AHe age-elevation show the 10% uncertainty considered for AHe ages. The catchment hypsometry frequency as well as glaciers and hillslope hypsometry frequencies are presented in (**b**). The bedrock AHe distribution considering 10% relative uncertainties and bedrock AFT distribution with > 30% relative uncertainties are shown in (**c**) and (**d**), respectively.

and 3.5) or accordingly (Sect. 3.4). This is similar to the proposed minimum number of ages to adequately resolve a component representing > 5% of the total detrital population (Vermeesch, 2004). We repeat this sampling process 10 000 times, storing the resulting SPDFs. From all of these detrital SPDFs, we compute the mean detrital SPDF and present the range of inferred detrital distributions.

# 3 Results

## 3.1 Ideal detrital age distributions

The bedrock age distributions for the two thermochronological systems, AHe and AFT (Fig. 4c–d), represent the ideal detrital age distributions if we sample the total catchment uniformly (allowing for the added noise). In the subsequent text, we refer to these ideal detrital age distributions as the bedrock distributions. The bedrock AHe SPDF shows a major peak at ~ 6.3 Ma, with the two youngest and smaller peaks at ~ 3.3 and ~ 4.3 Ma, as well as an older one at ~ 11 Ma. The shape of this SPDF mimics the hypsometric curve (black curve, Fig. 4b), which is expected with a uniform production of sediments and a linear age-elevation relationship. The AFT SPDF shows a single major peak

at ~ 16 Ma, and the form does not mimic the hypsometric curve. This partly reflects the form of the age-elevation relationship, in which the ages below 2000 m (Fig. 4a) are similar, limiting the ability to resolve a unique relationship of age and elevation. We also note that, due to the high uncertainty, the probability to have an AFT age equal to zero is non-null.

## 3.2 Transport of detrital particles

Here, we investigate the kinematics of particles in our reference model for the Tiedemann Glacier at the equilibrium state. To this end, we produce one particle in each cell of the model simultaneously (Fig. 5a) and let particles be transported away, potentially to the glacier front (i.e. the frontal moraine) where they would be deposited (Fig. 5b–d). In this case, assuming we have reached steady state in terms of transport, differences between the final detrital and the bedrock SPDFs can only be related to the process of particle transport.

We investigate the evolving contribution of particle source locations to the frontal moraine and to the inferred detrital age distributions (Fig. 6). Particles are sampled independently of their source origin (hillslope vs. glacial; Fig. 5a). According to our model, the minimum time required for a
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Kinematics of particle transport. Particles have been uniformly produced across the catchment and transported. (a) Sources of particles; red dots and blue dots are particles originating from the hillslope and glacier at the beginning of the model simulation, respectively. (b) Particle velocity at time 500 yr, (c) at time 1000 yr, and (d) at time 8500 yr. AHe ages in the frontal moraine (FM) are shown in (e); see also the inset in (d).</figcaption></figure>

particle originating from elevated parts of the catchment (i.e. far from the glacier front) to reach the frontal moraine is ~500 years (Fig. 6a, pink curve). Furthermore, the proportion of particle source locations in the frontal moraine becomes nearly constant after ~1500 years. After 8500 years, some particles are still close to their sources (Fig. 5d). We find that only 44% of the initial particles have reached the frontal moraine (Figs. 6b and S5a). The other 56% is trapped upstream, with a large portion of them resting in a lateral moraine (Figs. 5d and 6b) or on the side of small tributary glaciers having very low velocities (i.e. < 1 m yr⁻¹; Fig. 3c). These low velocities seem to result from a morphodynamic

feedback as the slope of tributary glaciers that carry the remaining particles is close to zero in the direction of sliding (Fig. S5b). This could, in turn, be related to the location of the ELA (Fig. 6b). As the number of particles within the frontal moraine does not increase significantly after 2000 years of simulation (i.e. reaching around 37% of the total number of initial particles; see Fig. S5a), we consider the detrital SPDF to have reached steady state by this time.

We observe a relative decrease in young ages (blue curves) in the detrital SPDFs (Fig. 6c–d) with time during the transient phase, reflecting the progressive arrival at the frontal moraine of older ages coming from upstream and elevated
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Transient evolution of particle sources and detrital thermochronological age distributions in the frontal moraine. (a) Cumulative contribution of source locations for particles and (b) their distribution with elevation after 8500 years of transport. Panels (c) and (d) show the transient evolution of the detrital thermochronological age distribution in the frontal moraine for AHe and AFT systems, respectively. BSPDF: bedrock SPDF.</figcaption></figure>

parts of the catchment (Fig. S6). The shape of detrital SPDFs reaches a steady state after around 2000 years of simulation.

The final form of the AHe detrital SPDF (i.e. red curve) shows a lower age peak at ~ 6.3 Ma and a larger proportion of ages between 8 and 12 Ma compared to the bedrock SPDF. The two youngest age peaks, ~ 3.3 and ~ 4.3 Ma, in the detrital SPDF are close to those of the bedrock SPDF. Thus, the results suggest that the ice transports and deposits mid-range ages, around 6 Ma (i.e. ~1500–2000 m), higher in the catchment.

For the AFT system (Fig. 6d), the magnitude of the age peak at ~ 16 Ma is also lower than that of the bedrock age distribution. The relative contribution of old ages (> 25 Ma) is a little higher in the detrital SPDF, but for ages younger than 10 Ma, the distributions are similar. The differences suggest that ice also transports and deposits sediments mostly originating from 1000–2000 m in this case. We note that the transient detrital SPDF at 1000 years (green curve, Fig. 6d) best matches the bedrock curve.

## 3.3 Detrital signature under uniform catchment erosion

Our model shows that more than half of the particles may not reach the frontal moraine. Consequently, we want to assess how the discrepancy between the detrital age distribution and the bedrock SPDF evolves if we consider the continuous production of particles (i.e. uniform erosion) across the catchment. To this end, we force the erosion rate to be 1 mm yr⁻¹ in each cell of the Tiedemann Glacier catchment (Fig. 7a). This is within the range of natural values, and it allows for the continuous production of particles while maintaining a reasonable simulation time (i.e. a particle is produced every 10 years in each grid cell). In this section, no distinction is made between the two sources in the sampling process when building detrital SPDFs. We stop the simulation after 2000 years, as the detrital SPDF is close to equilibrium (Fig. 6), and because the total number of particles increases rapidly, which increases the computational time. Here, we focus our analysis on the AHe system for clarity, but the detrital
SPDFs for the AFT system can be found in the Supplement (Fig. S7).

First, we observe a large variability of detrital SPDFs resulting from sampling 105 particles in the frontal moraine (Fig. 7c). However, the bedrock age distribution is always included in the range of inferred detrital SPDFs. Second, the mean detrital SPDF (red curve) represents the “true” detrital signal we expect to obtain by considering all the particles in the sampling site. This mean detrital SPDF differs from the previous model (Fig. 6c) in that the younger ages are overrepresented with two high peaks (at ~3.3 and ~4.3 Ma); also, there is an excess of older ages (8–11 Ma). However, mid-range ages of ~5 to ~8 Ma are still under-represented compared to the bedrock SPDF. Looking at the cumulative probability age distributions (CSPDFs), differences between the detrital CSPDF and the bedrock CSPDF are less clear, with the most obvious differences between 3 and 6 Ma.

The detrital AFT SPDFs (Fig. S7) show a similar tendency in under-representing mid-ages (~16 Ma) but much less pronounced, with a maximum difference with the bedrock SPDF at ~12 Ma and a smaller one at ~28 Ma. The AFT bedrock and detrital CSPDFs have a maximum difference at ~19 Ma but are similar overall. The results show that detrital SPDFs are shifted toward younger ages and deviate from what is expected for the case of uniform sampling (i.e. we expect the bedrock SPDF). We assign this behaviour to the short total transport distance for such ages. Furthermore, the SPDF is more informative on the effect of transport than the CSPDF that tends to smooth the signal. Finally, the quality of the results depends on the age uncertainty of the thermochronological system and on the age-elevation profile.

## 3.4 Detrital signature of glacier and hillslope sources

To understand the form of the mean detrital SPDF above, we now assess the respective influence of the two main sources of eroded particles: ice-free hillslopes and glaciers. We first consider the ice-free hillslope sources (Fig. 7a). Within the frontal moraine we sample only particles that originated from hillslopes. Thus, the detrital AHe ages range from 3.1 to 13.8 Ma, spanning ~96% of the total age range shown by the age-elevation profile (i.e. from 2.7 to 13.8 Ma). We calculate a new bedrock SPDF according to the hypometric distribution of hillslope sources and the AHe age-elevation profile (Fig. 4). The range of inferred detrital SPDFs is lower than the previous case (i.e. grey shaded area) but still remains important. The mean detrital SPDF shows a plateau-like signal from ~5 to ~9.5 Ma, suggesting that elevations from ~1400 to ~2500 m contribute in the same proportion to the particle budget in the frontal moraine. However, a comparison with the hillslope bedrock SPDF indicates that storage of sediment occurs for such elevations. Moreover, the hillslope mean detrital SPDF is shifted toward older ages (> 8 Ma) and excludes ages lower than 3.5 Ma as hillslope erosion does not operate at the lower elevations which have these young ages

(Fig. 4b). The hillslope mean detrital CSPDF shows a maximum difference with the hillslope bedrock CSPDF around 8 Ma (Fig. 8b). Therefore, constraining the hillslope sources for the bedrock SPDF allows us to better identify storage locations of sediment in comparison to the catchment bedrock SPDF (Fig. 8a–b).

The mean AFT detrital SPDF (Fig. S8a) shows a similar tendency of over-representing old ages (> 25 Ma) and under-representing mid-range ages (~16 Ma) compared to the hill-slope bedrock SPDF. However, the mean AFT detrital SPDF differs from the mean detrital AHe SPDF as it does not show a plateau-like form. Indeed, two age peak components occur at ~16 Ma and at ~40 Ma. However, the detrital cumulative distribution (Fig. S8b) is also shifted toward older ages, with a maximum difference at 24.5 Ma. The difference between the AHe and AFT age distributions is attributed to the different age-elevation relationships, with less spatial age resolution for the AFT system (Fig. 4a). Overall, the hillslope SPDFs reveal that the relevant source region ages tend to have older ages for both AHe and AFT, which is in accordance with the hillslope hypometric distribution (Fig. 4b).

We now consider particles originating from glacial sources (Fig. 7a), which we define to be all the ice-covered areas. The associated detrital AHe ages now range from 2.7 to 12.1 Ma, according to the hypometric distribution of the glacial source part of the age-elevation profile (Fig. 4), spanning the first ~85% of the total age range of the catchment. The maximum age reflects the maximum elevation reached by the ice (i.e. ~3420 m). As before, we calculate a new bedrock SPDF based on the glacial source distribution (Fig. 8). The detrital SPDFs, representing particles of glacial origin in the frontal moraine, show a maximum range of inferred SPDFs at ~3.3 Ma (Fig. 8c). We relate this range to the constant 10% relative age uncertainty applied, as the range of inferred detrital SPDFs tends to increase with younger ages, and younger ages have smaller absolute uncertainties. Compared to the glacial bedrock SPDF, the mean glacial detrital SPDF over-represents young ages (< 5 Ma) in contrast to the hillslope source, as we might expect. Next, we observe a relative under-representation of ages ~6.3 Ma and ages older than 8 Ma, consistent with a storage of sediments at mid-elevation. This is reflected in the cumulative detrital distribution (Fig. 8d), which is shifted towards younger ages and does not intersect the glacier bedrock cumulative distribution. The maximum difference between the two cumulative distributions is at 4.6 Ma.

The detrital AFT age distributions (Fig. S8) show the same tendency of over-representing young ages (age peak ~16.5 Ma) compared to the glacier bedrock distribution and under-representing old ages (> 28 Ma). The detrital CSPDF is also shifted toward younger ages; the maximum difference with the bedrock CSPDF is observed at ~28 Ma. We notice, however, smaller differences between the detrital and the bedrock distributions for the AFT than for the AHe system. Overall, the results show that older ages mainly reflect
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Spatial distribution of particles and detrital AHe age distributions within the frontal moraine. (a) Source locations of particles, with Er being the erosion rate. The spatial AHe age distributions of particles after 2000 years of transport (b). The detrital AHe age probability distributions (c) and their cumulative values (d). The grey shaded area represents the range of the inferred 10 000 detrital SPDFs from the sampling process. The black square in (b) shows the frontal moraine (FM) position.</figcaption></figure>

hillslope sources and younger ages the glacial sources, which is expected given the age-elevation profile and hypsometric curves of the two sources (Fig. 4b). However, they also show that the estimation of the bedrock distributions of the sources helps better constrain the role of ice transport and also better identify the predominant source in the sediment at the sampling site.

## 3.5 Detrital signature using a landscape evolution model with non-uniform erosion

To illustrate how the bedrock distributions of the sediment sources can help identify the dominant source signal at the sampling site, we now consider a model with non-uniform erosion. We use erosion laws presented in Sect. 2.1 with parameters listed in Table 1, which leads to the erosion pattern in Fig. 9a. Eroded particles are thus a mixture of hillslope- and glacier-origin sources. The parameters of the erosion laws were calibrated by a trial-error process to obtain a mean erosion rate over both source regions of $1 \text{ mm yr}^{-1}$ for a more consistent comparison with the constant-erosion-rate models. Obviously, variations occur locally about this mean value (e.g. maximum erosion rate is $\sim 31 \text{ mm yr}^{-1}$; Fig. 9a), and the calculated standard deviations for the mean erosion rates for the hillslope and glacier sources are $1 \pm 0.29$ and $1 \pm 2.9 \text{ mm yr}^{-1}$, respectively. Variations in ice erosion result mostly from the distribution of effective pressure that con-

trols the opening of cavities (Eqs. 3 and 6) and therefore from the length on which the basal shear stress is applied (Eq. 1 and Fig. S1). Consequently, there are some ice-covered areas that do not erode at all. The range of detrital ages for the glacier-origin particles spans 9 Myr from 2.7 to 11.7 Ma, whereas those for the hillslope-origin particles range from 3.5 to 13.8 Ma.

The mean detrital SPDF shows a large age peak at $\sim 5$ Ma and a smaller one at $\sim 10$ Ma. The glacial bedrock distribution predominates for ages < 5.5 Ma (Fig. 9c), suggesting a glacial source for the first age peak. This is confirmed by the high erosion rates located where the AHe ages are around 5 Ma (see Figs. 9a–b and S6a). The second age peaks (10 Ma) best match the hillslope bedrock distribution, which would suggest a dominant hillslope source contribution. This assumption is also confirmed by erosion rates of $\sim 1.80 \text{ mm yr}^{-1}$, with AHe ages at $\sim 10$ Ma (Figs. 9a and S6). The mean detrital SPDF also shows a low at $\sim 7–8$ Ma that we attribute to both the effect of storage at mid-elevation (Fig. 6b) and the erosion pattern, as the relevant elevations are represented by small tributary glaciers with low erosion rates ($< 0.07 \text{ mm yr}^{-1}$). A similar comparison can be done with the CSPDFs (Fig. 9d), although the curves are smoother. Thus, taking into account the bedrock age distributions of the different sources increases the potential to discriminate these sources in a detrital sample. However, we stress that the bedrock distributions are based on the assumption of spa-
<figure><img src="image_10.png" /></figure>

**Figure 8.** Detrital AHe age distributions of the frontal moraine and their cumulative probability for models considering only ice-free hillslope sources (**a–b**) and glacial sources (**c–d**). The grey shaded area represents the range of the inferred 10 000 detrital SPDFs from the sampling process. MSPDF: mean synoptic probability distribution.

<figure><img src="image_11.png" /></figure>

**Figure 9.** Non-uniform erosion model experiment with (**a**) the erosion rate, (**b**) spatial AHe age distribution of particles, (**c**) detrital AHe age probability distributions, and (**d**) their cumulative probability for the frontal moraine. The grey shaded area represents the range of the inferred 10 000 detrital SPDFs from the sampling process. The black square in (**b**) shows the frontal moraine (FM) location. MSPDF: mean synoptic probability distribution function.
tially uniform erosion across the catchment. Non-uniform erosion leads to different bedrock SPDF for the sources and also different expected sediment source distributions.

The mean detrital AFT SPDF (Fig. S9) is in agreement with a glacial source of particles for AFT ages of ~ 16 Ma, at which the erosion rates are highest (Fig. S9a), but this is much less pronounced than for the AHe system. Moreover, it is difficult to make a distinction between the two sources (hillslope vs. glacier) for ages > 30 Ma, as the mean detrital SPDFs closely match all the bedrock SPDFs. We relate this pattern to the large relative uncertainties (i.e. > 30%) we considered for the AFT ages that over-smooth the SPDF, especially for older ages. The mean detrital CSPDF shows no significant differences with the catchment bedrock CSPDF and suggests uniform erosion across the catchment, which is in contrast to the conclusion inferred for the AHe system, i.e. a major contribution of the glacier source. This difference, due to the different age-elevation relationships and age uncertainties between the two thermochronometers, reduces the ability to track the source of particles with AFT data in this example.

# 4 Discussion

## 4.1 Limitations of the particle transport model

It is important to note that the sediment particles in our model are passively transported and do not interact with the ice. For instance, the concentration of debris in the basal ice can influence the rheology of the ice and the friction of the bed, which both impact the ice flow (e.g. Hallet, 1981; Iverson et al., 2003; Cohen et al., 2005). Surface debris cover can also shield the ice from solar energy and, in turn, reduces the ablation rates and increases the local mass balance of the glacier (e.g. Östrem, 1959; Kayastha et al., 2000; Rowan et al., 2015). As iSOSIA does not perform a full 3D computation of ice deformation, the horizontal velocity of particles is parameterized as a simple polynomial function (i.e. Eq. 9) of the average ice velocity (i.e. for englacial particles). Sediments can move vertically only by the accumulation of new ice or by melting, whereas natural particles can also move due to internal ice deformation, including thrusting and folding (Hambrey et al., 1999; Hambrey and Lawson, 2000). Surface debris may also fall into crevasses and therefore be incorporated into the basal ice (Hambrey et al., 1999).

Additionally, our models omit debris comminution during transport. This process may have important implications for sampling strategies as block-sized sediments formed by e.g. quarrying potentially originate close to the sampling site (i.e. have a short travel distance) but are often not sampled. This could bias detrital SPDFs that are only representative of the sampled size fraction of sediments (mostly sand size) and thus sources located high in the catchment (i.e. larger travel distances). Lastly, sediment transport by rivers and subglacial drainage is neglected in the models presented here. However,

it has been shown that a large portion of basal debris can be evacuated by meltwater on a seasonal timescale (Collins, 1996; Kirkbride, 2002; Swift et al., 2005; Delaney et al., 2018). Not allowing for this factor potentially impacts the transfer time of our sediment particles to the glacier front, and we return to this point below. Additionally, water mixes sediments efficiently and may reduce the role of advective sediment transport in the ice. The role of sediment mixing by water has been investigated by Enkelmann et al. (2015), and they conclude that sampling the pro-glacial river (Ehlers et al., 2015) is similar to sampling through the ice-core terminal moraine, meaning that sediments in glacial outwash have a greater potential to be mixed by meltwater. We emphasize that our simple particle transport model is a starting point for studying the behaviour of sediment transport by ice, and integration of other processes of transport should be considered in future studies.

## 4.2 Effect of sediment transport by ice

Despite the limitations mentioned above, some first-order insights related to ice dynamics can be identified. Firstly, our ice transport model for the Tiedemann Glacier implies an equilibrium time in terms of the relative proportions of particle sources in the frontal moraine of the order of $10^3$ years (~1500 years). This timescale is of the order of the characteristic timescale developed by previous studies (Jóhannesson et al., 1989; Oerlemans, 2001; Roe and O'Neal, 2009; Herman et al., 2018) and mainly reflects the glacier dynamics. Incorporating transport of sediments by meltwater, this timescale would probably be reduced for glaciers with well-developed drainage systems, as a large portion of glacier basal debris seems to be evacuated by the subglacial drainage system (Collins, 1996; Kirkbride, 2002; Swift et al., 2005; Delaney et al., 2018). However, the types of debris forming frontal moraines vary from glacier to glacier. For instance, debris-covered glaciers are dominated by supraglacial debris that is mainly transported passively by ice (Benn and Evans, 2014). Terminal moraines of debris-covered glaciers are thus mainly built by a process of dumping debris from the ice surface and thus reflect glacier dynamics, especially if the glacier remains stable for a long period of time (Sharp, 1984; Lukas, 2005; Hambrey et al., 2008; Benn and Evans, 2014). Moreover, some sedimentological studies on terminal and lateral moraines have shown limited amounts of glaciofluvial-related facies (e.g. Winkler and Matthews, 2010; Bowman et al., 2018; Ewertowski and Tomczyk, 2020), as is the case for the Tiedemann Glacier (Menounos et al., 2013). Moreover, a recent study shows sediment transfer times relevant to building moraines of the order of 100–200 kyr (Cogez et al., 2018). Therefore, the timescale for building the frontal moraine would strongly depend on (i) the availability of subglacial debris vs. supraglacial debris and (ii) the glacier dynamics. We recommend that future studies focus on the timescales for build-
ing the terminal moraine, but we also highlight that, in some glacial areas, this timescale could impact the interpretation of detrital SPDFs from such glacial features, e.g. by reflecting a transient phase of sediment arrival (as shown in Fig. 6).

Secondly, zones of slow ice flow, such as in small tributary glaciers, may act as sediment traps as shown by particles remaining in elevated areas after 8500 years of simulation. These zones of slow ice motion occur mostly around the ELA and are associated with low values of local topographic slope. This is consistent with the suggestion that glacial erosion is more efficient around the ELA (e.g. Egholm et al., 2009; Brozović et al., 1997; Anderson et al., 2006; Steer et al., 2012), which leads to limited local relief and slope around the ELA and in turn to limited ice-sliding velocity. This also highlights morphodynamic feedbacks controlling, here by a trapping effect, the detrital distribution of thermochronological ages. Therefore, only ~44% of the total initial number of particles reach the frontal moraine after 8500 years of transport (i.e. the maximum time of our simulation). About 25% of the particles are stored in the lateral moraine (Fig. 5d); the remaining ~31% are trapped at higher elevations and therefore have residence times greater than 8500 years. The robustness of the results presented so far to different parameters has been tested by varying the glacier size and the hillslope diffusivity for particle transport. The results are presented in Sect. S2 of the Supplement. Overall, the conclusions are similar to those already discussed and show equilibrium times for the frontal moraine of the same order, with around 53%–60% of the total sediments stored higher in the catchment. To illustrate the effect of low velocities on transport times in more detail, we compute the average transfer time as a function of source location for particles that reach the frontal moraine (Fig. 10). The results show that the time required for a particle to reach the frontal moraine is not simply proportional to the distance from the source location. Indeed, some particles formed near the frontal moraine may take more than 3000 years to reach the glacier front due to velocities close to zero (Fig. 10). In contrast, some particles formed far from the frontal moraine have transfer times less than 1000 years due to averaged velocities greater than $15 \text{ m yr}^{-1}$ (Fig. 10b). The proximity of the source to ice streams with high velocities explains this pattern of transfer times (Fig. 3c). The average transfer times for particles formed along hillslopes or glaciers are $1825.4 \pm 1914.7$ and $1084.9 \pm 1014$ years, respectively, but show high variability. This difference is controlled (i) by the longer average distance of the hillslopes to the frontal moraine, $15.26 \pm 7.19$ km, compared to glaciers, $13.85 \pm 6.15$ km, and (ii) by spatial variation in ice flow velocities and storage of sediments in small tributary glaciers, which is more pronounced for the hillslope sources as shown in Fig. 10. Overall, sediment transfer times of our models are strongly influenced by the spatial distribution of small tributary glaciers, implying an important control of glacier sizes on the delivery of sediments to the main glacier transport system.

In the case of uniform erosion and a simple relationship between thermochronological age and elevation (Fig. 4a), we expect the detrital thermochronological signal associated with each source to mimic its associated hypsometric distribution. However, the resulting detrital SPDFs differ from those expected from the hypsometric distributions (Fig. 7). The lack of ages in the detrital SPDF corresponding to the peak observed ~1500–2000 m in the hypsometric curve may be explained by (i) the storage of a major portion of the sediments (i.e. 56%) outside the sampling site and (ii) by the patterns of transfer times for particles that reached the frontal moraine (Fig. 10). The generality of this conclusion is limited by our sediment transport model, as mentioned earlier, and the discrepancy observed between the detrital SPDFs and the expected distribution (i.e. bedrock SPDF) can be reduced for glaciers with a well-developed subglacial drainage system.

## 4.3 The role of detrital sources: glaciers or hillslopes

According to our model for the Tiedemann Glacier, each type of detrital source, i.e. hillslope or glacier, displays a different average hypsometry (Fig. 4b). Glaciers mostly represent the valley floor, while hillslopes mostly represent elevations greater than ~1800 m. Hillslope sources contribute to older ages in the detrital SPDFs, as expected from the age–elevation profiles and the range of elevation for hillslope sources (~700–4000 m). The glacial detrital SPDF (Fig. 8) mainly reflects the younger ages from the valley floor and a lack of older ages in the elevation range occupied by glaciers (~530–3420 m).

When estimating glacial erosion rates from sediment flux measurements it is important to distinguish glacier-origin debris from supraglacial debris. Previous studies used cosmogenic nuclides concentrations (e.g. Guillon et al., 2015) or U–Pb ages (e.g. Godon et al., 2013) on sediments to discriminate between sources. Here, we tested a simple approach by characterizing the bedrock age distribution of different sources and compared them with the detrital age distribution. The model considering non-uniform erosion gave promising results for identifying the contribution of different sediment sources to the frontal moraine (Fig. 9).

Despite identical averaged erosion rates for the two sources (i.e. $1 \text{ mm yr}^{-1}$), the resulting shape of the mean detrital SPDF is the result of higher local glacial erosion rates compared to the more diffusive erosion on hillslopes. Therefore, the glacial source locally produces a high quantity of sediment particles with similar thermochronological ages, which are ultimately transported to the frontal moraine. This explains the high peak observed at ~5 Ma (Fig. 9) and the lower peak at ~10 Ma corresponding to hillslope sources. Furthermore, proximity to the depositional site of the sources with high glacier velocities also contributes to the magnitude of the age peak components observed in the detrital SPDF (Fig. 10). Therefore, detrital SPDFs of glacial features such as the frontal moraine are likely to over-represent
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Characteristics of particle transfer to the frontal moraine: (a) the time needed for each particle to reach the frontal moraine and (b) their averaged speed. Black dots are particles that did not reach the frontal moraine (FM).</figcaption></figure>

glacial sources compared to hillslope sources if driven by locally high erosion rates. Finally, a large proportion of sediment particles produced on hillslopes originate from around 2800 m as illustrated by the older age peak component at ~10 Ma (Fig. 9c). However, in our models the hillslope sources are more subject to storage by tributary glaciers, which reduces their contribution to the detrital signal in the frontal moraine.

## 4.4 Implications for detrital sampling strategies

A sampling strategy equivalent to the one considered in this study (i.e. randomly sampling the entire frontal moraine) has the potential to capture more bedrock age components, as proposed in previous studies (e.g. Enkelmann et al. and Ehlers, 2015). To illustrate this effect, we now consider different sampling strategies. We perform an additional experiment with uniform erosion, wherein sampling occurs in four different regions of the frontal moraine (Fig. 11). The process of sampling is the same as for previous models; i.e. we randomly collect 105 particles within each region, produce an SPDF, and repeat this process 10 000 times to infer a mean detrital SPDF. The resulting mean AHe detrital SPDFs and CSPDFs (Fig. 11) show significant variability. Sampling region 1 mostly captures young ages (< 6 Ma) and therefore the glacier source, while hillslope sources with older ages

(> 6 Ma) are under-represented. Sampling region 2, located closer to centre of the moraine, includes an older age component (> 6 Ma), which leads to a better fit with the mean detrital SPDF obtained by sampling the entire frontal moraine. Statistical tests performed on the CSPDFs confirm this similarity (see Table A1). However, the different sensitivity of these statistical tests complicates the interpretation of these results. Finally, regions 3 and 4 show similar distributions, with an over-representation of young (< 4 Ma) and old ages (> 7 Ma), as well as a gap in mid-range ages (between 5 and 7 Ma) representing intermediate elevations (1500–2000 m). Old ages (red dots in Fig. 11a–b) are mostly transported in the central part of the main glacier and deposited in the central part of the frontal moraine (regions 2, 3, and 4). An exception occurs with regions 3 and 4 due to deviation of ice streamlines (Fig. 11b). This deviation is also responsible for significant deposition of sediments in the lateral moraine upstream. The trend of decreasing young age components in the detrital SPDF across the frontal moraine is seen for the AFT system, although the differences are less pronounced. However, this trend seems to be supported by a comparison with true ice-cored terminal moraine AFT data from Enkelmann and Ehlers (2015) (see Fig. B1). We observe that, from left to right, the detrital age distributions within the frontal moraine incorporate older ages. However, we are aware that the method used to build SPDFs with a Gaussian kernel
<figure><img src="image_13.png" /></figure>

**Figure 11.** Detrital AHe age distributions of the five regions seen in (b) (black rectangles) from the experiment considering a uniform source of particles. Spatial distribution of particles with their associated AHe ages (a), with a zoom-in of the frontal moraine area in (b). The density plots and their cumulative distributions are shown in (c) and (d). The dashed black line is the mean detrital SPDF for the frontal moraine (FM).

(Eq. 11) tends to break down with high (> 30%) relative uncer- tainties (Brandon, 1996), as for the presented AFT data. We use this method for simplicity and to facilitate compar- ison with the synthetic AFT SPDFs. However, the original data also show the same tendency of age components get- ting older across the terminal moraine, and statistical tests applied to original data (Enkelmann and Ehlers, 2015) sup- port the high variability of the detrital AFT age distributions across the terminal moraine regions (see Table A1). Overall, local sampling within the frontal moraine or upstream along a transverse transect leads to a higher variability in the in- ferred detrital age distributions. Even if some regions may show better agreement with the bedrock SPDF (region 2), randomly sampling small patches of sediments through the moraine captures most of the age components on average and seems to be a better strategy overall.

## 4.5 The effect of age uncertainties and age-elevation profiles: comparing AHe and AFT

For all of our models, we concluded that detrital age dis- tributions resulting from the AFT ages were less informa- tive than those from the AHe ages. The differences between the two systems occur for two main reasons. First, the age- elevation profiles differ. For the AFT profile, the youngest ages are not at the lowest elevation but occur at ~ 2000 m of elevation. The slope of the AFT age-elevation relationship is

almost vertical (and actually negative) for elevations lower than 2000 m. This low age-elevation gradient leads to a re- duction in source identifiability; i.e. similar ages can come from a large range of elevations. Secondly, the high relative age uncertainties (i.e. > 30%) in the AFT data smooth the SPDF and decrease the resolvability of age components in the age distribution. Consequently, these two characteristics can make the AFT system perhaps less useful for tracking erosion patterns, as seen for the case of non-uniform erosion. However, in some cases the AFT distribution can still capture first-order behaviour of sources if combined with bedrock age and age-elevation distributions, depending on the distri- bution of erosion across the catchment.

## 5 Conclusions

In this study, we have presented a numerical approach to investigate the effect of sediment sources (hillslopes vs. glaciers), ice transport, and deposition on the distribution of thermochronological ages found in a frontal moraine. We ap- plied this approach to a glaciated catchment which presents simple morphological characteristics: the Tiedemann Glacier (British Columbia, Canada) in steady state. Firstly, small tributary glaciers with very low velocities may act as traps for sediments and delay particle transfer to the glacier front. These low velocities may result from a morphodynamic feed- back between the location of the ELA and the bed slope in
the direction of sliding. Secondly, the transfer times of sedi-
ments are influenced by the proximity of their sources to the
ice streams showing high velocities. Indeed, sediments lo-
cated in elevated areas may experience lower transfer times
than sediments produced closer to the sampling site. Thirdly,
horizontal separation of ice flow lines can produce lateral
moraines that may store a significant amount of sediment
produced upstream. This implies that frontal moraines may
include sediments that contain thermochronological signa-
tures from limited parts of the total catchment. To address
this problem, lateral moraines of the same age deposition as
the frontal moraine can also be targeted for complementary
sampling, therefore incorporating more age components.

Sediment transport by ice can lead to differences in detri-
tal age distributions compared to the bedrock age distribu-
tion, for instance by undersampling mid-altitude age compo-
nents from tributary glaciers showing low sliding velocities
(< 1 m yr⁻¹). This could lead to misinterpretation of regional
erosion patterns. Moreover, strategies considering local sam-
pling of sediment in the frontal moraine show variable detri-
tal age distributions that predominantly reflect the variabil-
ity of local sources upstream. In principle, this may allow
us to directly associate particle sinks, e.g. moraines, with
their sources. In contrast, randomly sampling through the
frontal moraine potentially captures more age components,
providing a more representative picture of the whole catch-
ment. Therefore, we suggest the sampling strategy should
be designed according to the question being addressed. Fur-
thermore, we systematically compared two thermochrono-
logical systems, AHe and AFT, with different but coherent
age-elevation profiles and different relative age uncertain-
ties. While the first factor plays a role in the ability to track
sediment sources, the second factor impacts the precision of
SPDFs. However, characterization of the spatial variability
of source contributions between different sampling sites may
still be possible depending on the distribution of erosion.

Overall, our numerical approach offers novel insights on
the application of detrital thermochronology to glaciated
catchments and on the role of long-term exhumation, mod-
ern erosion, and past sediment transport. However, as we
stated earlier, this study considers simple laws for ice motion
and sediment entrainment and neglects sediment transport by
the glaciofluvial system. Clearly, directions for future studies
include the role of subglacial hydrology and the plastic de-
formation of ice in sediment transport as well as the role of
meltwater in the building of terminal moraine. Finally, we
have only considered one single alpine catchment, and our
results may not be applicable to catchments showing more
complexity (e.g. tributary glacier valleys).
## Appendix A

**Table A1.** Statistical test results from Kuiper and Kolmogorov-Smirnov (KS) tests, with associated *p* values for the modelled detrital CSPDFs shown in Sect. 4.4. Each modelled detrital CSPDF (regions 1–4 and frontal moraine) is tested against the modelled catchment bedrock CSPDF. “Frontal moraine” corresponds to the mean detrital CSPDF of the entire moraine. The variability of detrital SPDFs of the original data from Enkelmann and Ehlers (2015), *p* values for S10–S14, has been tested by comparing the original ice-cored detrital SPDFs against the detrital age distribution of the glacial outwash sample (9TETG15) presented in Ehlers et al. (2015). Highlighted in black are contradictions or *p* values close to the alpha level (0.05) between the two statistical tests about the similarity between the corresponding detrital CSPDF and the catchment bedrock CSPDF (or glacial outwash in case of S10–S14). 1: the two distributions are different, 0: the two distributions are similar.

<table><thead><tr><th rowspan="2"></th><th colspan="2">AHe</th><th colspan="2">AFT</th></tr><tr><th>Kuiper test / p value</th><th>KS test / p value</th><th>Kuiper test / p value</th><th>KS test / p value</th></tr></thead><tbody><tr><td>Region 1</td><td>1 / 1.16 × 10<sup>-5</sup></td><td>1 / 4.10 × 10<sup>-7</sup></td><td>1 / 1.37 × 10<sup>-6</sup></td><td>1 / 3.25 × 10<sup>-8</sup></td></tr><tr><td>Region 2</td><td><strong>0 / 0.057</strong></td><td><strong>0 / 0.26</strong></td><td><strong>0 / 0.09</strong></td><td><strong>1 / 0.01</strong></td></tr><tr><td>Region 3</td><td>1 / 1.63 × 10<sup>-9</sup></td><td>1 / 4.80 × 10<sup>-8</sup></td><td>1 / 6.42 × 10<sup>-22</sup></td><td>1 / 3.26 × 10<sup>-24</sup></td></tr><tr><td>Region 4</td><td>1 / 4.46 × 10<sup>-11</sup></td><td>1 / 4.68 × 10<sup>-9</sup></td><td>1 / 5.78 × 10<sup>-20</sup></td><td>1 / 3.38 × 10<sup>-22</sup></td></tr><tr><td>Frontal moraine</td><td><strong>0 / 0.50</strong></td><td><strong>0 / 0.55</strong></td><td>0 / 0.87</td><td>0 / 0.31</td></tr><tr><td>S10</td><td>—</td><td>—</td><td>1 / 8.83 × 10<sup>-16</sup></td><td>1 / 1.30 × 10<sup>-17</sup></td></tr><tr><td>S11</td><td>—</td><td>—</td><td>1 / 4.99 × 10<sup>-9</sup></td><td>1 / 4.03 × 10<sup>-7</sup></td></tr><tr><td>S12</td><td>—</td><td>—</td><td><strong>1 / 0.03</strong></td><td><strong>0 / 0.30</strong></td></tr><tr><td>S13</td><td>—</td><td>—</td><td>0 / 0.25</td><td>0 / 0.38</td></tr><tr><td>S14</td><td>—</td><td>—</td><td>0 / 0.52</td><td>0 / 0.10</td></tr></tbody></table>

## Appendix B

<figure><img src="image_14.png" /></figure>

**Figure B1.** Detrital AFT age distributions of the four regions seen in (b) (black rectangles) from the experiment considering a uniform source of particles (a). The density plots of modelled detrital AFT ages are shown in (c). The dashed black line is the mean detrital SPDF for the frontal moraine (FM). The ice-cored detrital SPDFs from Enkelmann and Ehlers (2015) are shown in (d). Each sample (S10–S14) is located in (b) by the white stars, and the age distributions have been built using the method explained in Sect. 2.4. The bedrock SPDF in (d) results from the bedrock single-grain age distribution presented in Enkelmann and Ehlers (2015).
**Code availability.** The iSOSIA version (iSOSIA_3.4.7b) and the external routine used to compute the detrital age distributions are publicly available here: https://github.com/davidlundbek/iSOSIA. (David Lundbeck Egholm and Maxime Bernard, 2020.) (https://help.github.com/en/github/creating-cloning-and-archiving-repositories/referencing-and-citing-content).

**Supplement.** The supplement related to this article is available online at: https://doi.org/10.5194/esurf-8-931-2020-supplement.

**Author contributions.** MB developed the model code for the computation of thermochronological ages and SPDFs, designed the experiments, and prepared the paper with contributions from PS and KG. DLE provided the model iSOSIA with the associated knowledge to allow its use. He also contributed to the final version of the original draft.

**Competing interests.** The authors declare that they have no conflict of interest.

**Acknowledgements.** We thank Todd Ehlers and the anonymous reviewer for their reviews that helped to improve the quality and clarity of this study. The two handling editors, Jean Braun and Andreas Lang, are also acknowledged for their guidance. We also thank Eva Enkelmann for having shared the AFT data from the Tiedemann Glacier, which have been used to produce synthetic AFT ages. We give special thanks to Peter Van der Beek, Stephane Bonnet, Benjamin Guillaume, and Pierre Valla for their advice and guidance that contributed to this study. Finally, we thank the Institut Français du Danemark for their financial support.

**Financial support.** This research has been supported by the Institut Français du Danemark (IFD).

**Review statement.** This paper was edited by Jean Braun and reviewed by Todd A. Ehlers and one anonymous referee.

## References

Alley, R. B., Cuffey, K. M., Evenson, E. B., Strasser, J. C., Lawson, D. E., and Larson, G. J.: How glaciers entrain and transport basal sediment: physical constraints, Quaternary Sci. Rev., 16, 1017–1038, https://doi.org/10.1016/S0277-3791(97)00034-6, 1997.

Anderson, R. S.: Evolution of lumpy glacial landscapes, Geology, 42, 679–682, https://doi.org/10.1130/G35537.1, 2014. Anderson, R. S., Molnar, P., and Kessler, M. A.: Features of glacial valley profiles simply explained, J. Geophys. Res.-Earth, 111, 1–14, https://doi.org/10.1029/2005JF000344, 2006.

Andrews, D. J. and Bucknam, R. C.: Fitting degradation of shoreline scarps by a nonlinear diffusion model, J. Geophys. Res.-Sol. Ea., 92, 12857–12867, 1987.

Avdeev, B., Niemi, N. A., and Clark, M. K.: Doing more with less: Bayesian estimation of erosion models with detrital thermochronometric data, Earth Planet. Sc. Lett., 305, 385–395, https://doi.org/10.1016/j.epsl.2011.03.020, 2011.

Beaumont C., Fullsack P., Hamilton J., McClay K. R. (Eds): Erosional control of active compressional orogens, Thrust Tectonics, Springer, Dordrecht, 1–18, https://doi.org/10.1007/978-94-011-3066-0_1, 1992.

Benn, D. and Evans, D. J.: Glaciers and glaciation, second edition, Hodder education, Routledge, London, England, pp. 816, https://doi.org/10.4324/9780203785010, 2014.

Bernard, T., Steer, P., Gallagher, K., Szulc, A., Whitham, A., and Johnson, C.: Evidence for Eocene-Oligocene glaciation in the landscape of the East Greenland margin, Geology, 44, 895–898, https://doi.org/10.1130/G38248.1, 2016.

Boulton, G. S., Coates D. R. (Eds.): Glacial Geomorphology: Processes and Patterns of Glacial Erosion, Springer, Dordrecht, 41–87, https://doi.org/10.1007/978-94-011-6491-7_2, 1982.

Bowman, D., Eyles, C. H., Narro-Pérez, R., and Vargas, R.: Sedimentology and Structure of the Lake Palcacocha Laterofrontal Moraine Complex in the Cordillera Blanca, Peru, Revista de Glaciares y Ecosistemas de Montaña, 5, 16–16, https://doi.org/10.36580/rgem.i5.27-42, 2018.

Brædstrup, C. F., Damsgaard, A., and Egholm, D. L.: Ice-sheet modelling accelerated by graphics cards, Comput. Geosci., 72, 210–220, https://doi.org/10.1016/j.cageo.2014.07.019, 2014.

Brædstrup, C. F., Egholm, D. L., Ugelvig, S. V., and Pedersen, V. K.: Basal shear stress under alpine glaciers: insights from experiments using the iSOSIA and Elmer/Ice models, Earth Surf. Dynam., 4, 159–174, https://doi.org/10.5194/esurf-4-159-2016, 2016.

Brandon, M. T.: Probability density plot for fission-track grain-age samples. Radiat. Meas., 26, 663–676, https://doi.org/10.1016/S1350-4487(97)82880-6, 1996.

Brewer, I. D., Burbank, D. W., and Hodges, K. V.: Modelling detrital cooling-age populations: insights from two Himalayan catchments, Basin Res., 15, 305–320, https://doi.org/10.1046/j.1365-2117.2003.00211.x, 2003.

Brozović, N., Burbank, D. W., and Meigs, A. J.: Climatic limits on landscape development in the northwestern Himalaya, Science, 276, 571–574, https://doi.org/10.1126/science.276.5312.571, 1997.

Champagnac, J. D., Valla, P. G., and Herman, F.: Late-Cenozoic relief evolution under evolving climate: A review, Tectonophysics, 614, 44–65, https://doi.org/10.1016/j.tecto.2013.11.037, 2014.

Clarke, G. K.: A short history of scientific investigations on glaciers, J. Glaciol., 33, 4–24, https://doi.org/10.3189/S0022143000215785, 1987.

Cogez, A., Herman, F., Pelt, É., Reuschlé, T., Morvan, G., Darvill, C. M., Norton, K. P., Christl, M., Märki, L., and Chabaux, F.: U-Th and ¹⁰Be constraints on sediment recycling in proglacial settings, Lago Buenos Aires, Patagonia, Earth Surf. Dynam., 6, 121–140, https://doi.org/10.5194/esurf-6-121-2018, 2018.

Cohen, D., Iverson, N. R., Hooyer, T. S., Fischer, U. H., Jackson, M., and Moore, P. L.: Debris-bed friction of
hard-bedded glaciers. J. Geophys. Res-Earth, 110, F02007, https://doi.org/10.1029/2004JF000228, 2005.

Cohen, D., Hooyer, T. S., Iverson, N. R., Thomason, J. F., and Jackson, M.: Role of transient water pressure in quarrying: A subglacial experiment using acoustic emissions, J. Geophys. Res.-Earth, 111, F03006, https://doi.org/10.1029/2005JF000439, 2006.

Collins, D. N.: Sediment transport from glacierized basins in the In Erosion and Sediment Yield: Global and Regional Perspectives: Proceedings of an International Symposium Held at Exeter, UK, from 15 to 19 July 1996, (No. 236, p. 85), IAHS, https://doi.org/10.3389/feart.2018.00175, 1996.

Delaney, I., Bauder, A., Werder, M. A., and Farinotti, D.: Regional and annual variability in subglacial sediment transport by water for two glaciers in the Swiss Alps, Front. Earth Sci., 6, 175, https://doi.org/10.3389/feart.2018.00175, 2018.

Egholm, D. and Bernard, M.: iSOSIA_3.4.7b, Zenodo, https://doi.org/10.5281/zenodo.3875297, 2020.

Egholm, D. L., Knudsen, M. F., Clark, C. D., and Lesemann, J. E.: Modeling the flow of glaciers in steep terrains: The integrated second-order shallow ice approximation (iSOSIA), J. Geophys. Res.-Earth, 116, F02012, https://doi.org/10.1029/2010JF001900, 2011.

Egholm, D. L., Nielsen, S., Pedersen, V. K., and Lesemann, J. E.: Glacial effects limiting mountain height, Nature, 460, 884–887, https://doi.org/10.1038/nature08263, 2009.

Egholm, D. L., Pedersen, V. K., Knudsen, M. F., and Larsen, N. K.: Coupling the flow of ice, water, and sediment in a glacial landscape evolution model, Geomorphology, 141, 47–66, https://doi.org/10.1016/j.geomorph.2011.12.019, 2012a.

Egholm, D. L., Pedersen, V. K., Knudsen, M. F., and Larsen, N. K.: On the importance of higher order ice dynamics for glacial landscape evolution, Geomorphology, 141, 67–80, https://doi.org/10.1016/j.geomorph.2011.12.020, 2012b.

Ehlers, T. A.: Crustal thermal processes and the interpretation of thermochronometer data, Rev. Mineral. Geochem., 58, 315–350, https://doi.org/10.2138/rmg.2005.58.12, 2005.

Ehlers, T. A., Szameitat, A., Enkelmann, E., Yanites, B. J., and Woodsworth, G. J.: Identifying spatial variations in glacial catchment erosion with detrital thermochronology, J. Geophys. Res.-Earth, 120, 1023–1039, https://doi.org/10.1002/2014JF003432, 2015.

Enkelmann, E. and Ehlers, T. A.: Evaluation of detrital thermochronology for quantification of glacial catchment denudation and sediment mixing, Chem. Geol., 411, 299–309, https://doi.org/10.1016/j.chemgeo.2015.07.018, 2015.

Enkelmann, E., Zeitler, P. K., Pavlis, T. L., Garver, J. I., and Ridgway, K. D.: Intense localized rock uplift and erosion in the St. Elias orogen of Alaska, Nat. Geosci., 2, 360–363, https://doi.org/10.1038/ngeo502, 2009.

Ewertowski, M. W. and Tomczyk, A. M.: Reactivation of temporarily stabilized ice-cored moraines in front of polythermal glaciers: Gravitational mass movements as the most important geomorphological agents for the redistribution of sediments (a case study from Ebbabreen and Ragnarbreen, Svalbard), Geomorphology, 350, 106952, https://doi.org/10.1016/j.geomorph.2019.106952, 2020.

Falkowski, S., Enkelmann, E., Drost, K., Pfänder, J. A., Stübner, K., and Ehlers, T. A.: Cooling history of the St. Elias syn-

taxis, southeast Alaska, revealed by geochronology and thermochronology of cobble-sized glacial detritus, Tectonics, 35, 447–468, https://doi.org/10.1002/2015TC004086, 2016.

Farinotti, D., Brinkerhoff, D. J., Clarke, G. K. C., Fürst, J. J., Frey, H., Gantayat, P., Gillet-Chaulet, F., Girard, C., Huss, M., Leclercq, P. W., Linsbauer, A., Machguth, H., Martin, C., Maussion, F., Morlighem, M., Mosbeux, C., Pandit, A., Portmann, A., Rabatel, A., Ramsankaran, R., Reerink, T. J., Sanchez, O., Stentoft, P. A., Singh Kumari, S., van Pelt, W. J. J., Anderson, B., Benham, T., Binder, D., Dowdeswell, J. A., Fischer, A., Helfricht, K., Kutuzov, S., Lavrentiev, I., McNabb, R., Gudmundsson, G. H., Li, H., and Andreassen, L. M.: How accurate are estimates of glacier ice thickness? Results from ITMIX, the Ice Thickness Models Intercomparison eXperiment, The Cryosphere, 11, 949–970, https://doi.org/10.5194/tc-11-949-2017, 2017.

Farinotti, D., Huss, M., Fürst, J. J., Landmann, J., Machguth, H., Maussion, F., and Pandit, A.: A consensus estimates for the ice thickness distribution of all glaciers on Earth, Nat. Geosci., 12, 168–173, https://doi.org/10.1038/s41561-019-0300-3, 2019.

Farley, K. A. and Stockli, D. F.: (U–Th)/He dating of phosphates: Apatite, monazite, and xenotime, Rev. Mineral Geochem., 48, 559–577, https://doi.org/10.2138/rmg.2002.48.15, 2002.

Fitzgerald, P. G. and Stump, E.: Cretaceous and Cenozoic episodic denudation of the Transantarctic Mountains, Antarctica: New constraints from apatite fission track thermochronology in the Scott Glacier region, J. Geophys. Res.-Sol Ea., 102, 7747–7765, https://doi.org/10.1029/96JB03898, 1997.

Galbraith, R. F.: Statistics for fission track analysis, 1st Edition, Chapman and Hall/CRC, New York, pp. 240, https://doi.org/10.1201/9781420034929, 2005.

Gallagher, K.: Evolving temperature histories from apatite fission-track data, Earth Planet. Sc. Lett., 136, 421–435, https://doi.org/10.1016/0012-821X(95)00197-K, 1995.

Gallagher, K. and Parra, M.: A new approach to thermal history modelling with detrital low temperature thermochronological data, Earth Planet. Sc. Lett., 529, 115872, https://doi.org/10.1016/j.epsl.2019.115872, 2020.

Glen, J. W.: Experiments on the deformation of ice, J. Glaciol., 2, 111–114, https://doi.org/10.3189/S002214300034067, 1952.

Glotzbach, C., Busschers, F. S., and Winsemann, J.: Detrital thermochronology of Rhine, Elbe and Meuse river sediment (Central Europe): implications for provenance, erosion and mineral fertility, Int. J. of Earth Sci., 107, 459–479, https://doi.org/10.1007/s00531-017-1502-9, 2018.

Godon, C., Mugnier, J. L., Fallourd, R., Paquette, J. L., Pohl, A., and Buoncristiani, J. F.: The Bossons glacier protects Europe's summit from erosion, Earth Planet. Sc. Lett., 375, 135–147, https://doi.org/10.1016/j.epsl.2013.05.018, 2013.

Goodsell, B., Hambrey, M. J., and Glasser, N. F.: Debris transport in a temperate valley glacier: Haut Glacier d'Arolla, Valais, Switzerland. J. Glaciol., 51, 139–146, https://doi.org/10.3189/17275605781829647, 2005.

Grabowski, D. M., Enkelmann, E., and Ehlers, T. A.: Spatial extent of rapid denudation in the glaciated St. Elias syntaxis region, SE Alaska, J. Geophys. Res.-Earth, 118, 1921–1938, https://doi.org/10.1002/jgrf.20136, 2013.

Guillon, H., Mugnier, J. L., Buoncristiani, J. F., Carcaillet, J., Godon, C., Prud'Homme, C., Van der Beek, P., and Vas-
sallo, R.: Improved discrimination of subglacial and periglacial erosion using ¹⁰Be concentration measurements in subglacial and supraglacial sediment load of the Bossons glacier (Mont Blanc massif, France), Earth Surf. Proc. Land, 40, 1202–1215, https://doi.org/10.1002/esp.3713, 2015.

Hallet, B.: Subglacial regulation water film, J. Glaciol., 23, 321–334, https://doi.org/10.3189/S0022143000029932, 1979.

Hallet, B.: Glacial abrasion and sliding: their dependence on the debris concentration in basal ice, Ann. Glaciol., 2, 23–28, https://doi.org/10.3189/172756481794352487, 1981.

Hallet, B., Hunter, L., and Bogen, J.: Rates of erosion and sediment evacuation by glaciers: A review of field data and their implications, Global Planet. Change, 12, 213–235, 1996.

Hambrey, M. J., Bennett, M. R., Dowdeswell, J. A., Glasser, N. F., and Huddart, D.: Debris entrainment and transfer in polythermal valley glaciers, J. Glaciol., 45, 69–86, https://doi.org/10.3189/S0022143000003051, 1999.

Hambrey, M. J., and Lawson, W.: Structural styles and deformation fields in glaciers: a review, Geol. Soc., London, Special Publications, 176, 59–83, https://doi.org/10.1144/GSL.SP.2000.176.01.06, 2000.

Hambrey, M. J., Quincey, D. J., Glasser, N. F., Reynolds, J. M., Richardson, S. J., and Clemmens, S.: Sedimentological, geomorphological and dynamic context of debris-mantled glaciers, Mount Everest (Sagarmatha) region, Nepal, Quaternary Sci. Rev., 27, 2361–2389, https://doi.org/10.1016/j.quascirev.2008.08.010, 2008.

Herman, F., Beaud, F., Champagnac, J. D., Lemieux, J. M., and Sternai, P.: Glacial hydrology and erosion patterns: A mechanism for carving glacial valleys, Earth Planet. Sci. Lett., 310, 498–508, https://doi.org/10.1016/j.epsl.2011.08.022, 2011.

Herman, F., Braun, J., Deal, E., and Prasicek, G.: The response time of glacial erosion, J. Geophys. Res.-Earth, 123, 801–817, https://doi.org/10.1002/2017JF004586, 2018.

Herman, F., Seward, D., Valla, P. G., Carter, A., Kohn, B., Willett, S. D., and Ehlers, T. A.: Worldwide acceleration of mountain erosion under a cooling climate, Nature, 504, 423–426, https://doi.org/10.1038/nature12877, 2013.

Hindmarsh, R. C.: The role of membrane-like stresses in determining the stability and sensitivity of the Antarctic ice sheets: back pressure and grounding line motion, Philos. T. Roy. Soc. A, 364, 1733–1767, https://doi.org/10.1098/rsta.2006.1797, 2006.

Hurford, A. J. and Green, P. F.: The zeta age calibration of fission-track dating, Chem. Geol., 41, 285–317, https://doi.org/10.1016/S0009-2541(83)80026-6, 1983.

Iverson, N. R.: A theory of glacial quarrying for landscape evolution models, Geology, 40, 679–682, https://doi.org/10.1130/G33079.1, 2012.

Iverson, N. R., Cohen, D., Hooyer, T. S., Fischer, U. H., Jackson, M., Moore, P. L., and Kohler, J.: Effects of basal debris on glacier flow, Science, 301, 81–84, https://doi.org/10.1126/science.1083086, 2003.

Jarvis, A., Reuter, R. I., Nelson, A., and Guevara, E.: Hole-filled seamless SRTM data V4, International Centre for Tropical Agriculture (CIAT), available at: http://srtm.csi.cgiar.org, (last access: 1 December 2019), 2008.

Jóhannesson, T., Raymond, C. F., and Waddington, E. D.: A simple method for determining the response time of glaciers, in: Oerlemans J. (Eds.) Glacier fluctuations and climatic change,

pp. 343–352, Springer, Dordrecht, https://doi.org/10.1007/978-94-015-7823-3_22, 1989.

Kayastha, R. B., Takeuchi, Y., Nakawo, M., and Ageta, Y.: Practical prediction of ice melting beneath various thickness of debris cover on Khumbu Glacier, Nepal, using a positive degree-day factor, IAHS-AISH P. 264, 71–81, 2000.

Kessler, M. A., Anderson, R. S., and Briner J. P.: Fjord insertion into continental margins driven by topographic steering of ice, Nat. Geosci., 1, 365–369, https://doi.org/10.1038/ngeo201, 2008.

Kirkbride, M. P.: Processes of glacial transportation, in: Modern and Past Glacial Environments, edited by: Menzies J., Butterworth-Heinemann, Canada, Butterworth-Heinemann, 147–169, https://doi.org/10.1016/B978-075064226-2/50009-X, 2002.

Koppes, M., Hallet, B., Rignot, E., Mouginot, J., Wellner, J. S., and Boldt, K.: Observed latitudinal variations in erosion as a function of glacier dynamics, Nature, 526, 100–103, https://doi.org/10.1038/nature15385, 2015.

Koppes, M. N. and Montgomery, D. R.: The relative efficacy of fluidial and glacial erosion over modern to orogenic timescales, Nat. Geosci., 2, 644–647, https://doi.org/10.1038/ngeo616, 2009.

Lukas, S.: A test of the englacial thrusting hypothesis of ‘hummocky’ moraine formation: case studies from the northwest Highlands, Scotland, Boreas, 34, 287–307, https://doi.org/10.1111/j.1502-3885.2005.tb01102.x, 2005.

Matsuoka, N. and Murton, J.: Frost weathering: recent advances and future directions, Permafrost Periglac., 19, 195–210, https://doi.org/10.1002/ppp.620, 2008.

MacGregor, K. C., Anderson, R. S., Anderson, S. P., and Waddington, E. D.: Numerical simulations of longitudinal profile evolution of glacial valleys, Geology, 28, 1031–1034, https://doi.org/10.1130/0091-7613(2000)28<1031:NSOGLP>2.0.CO;2, 2000.

MacGregor, K. R., Anderson, R. S., and Waddington, E. D.: Numerical modelling of glacial erosion and headwall processes in alpine valleys, Geomorphology, 103, 189–204, https://doi.org/10.1016/j geomorph.2008.04.022, 2009.

Menounos, B., Clague, J. J., Clarke, G. K., Marcott, S. A., Osborn, G., Clark, P. U., Tennant, C., and Novak, A. M.: Did rock avalanche deposits modulate the late Holocene advance of Tiedemann Glacier, southern Coast Mountains, British Columbia, Canada?, Earth Planet. Sc. Lett., 384, 154–164, https://doi.org/10.1016/j.epsl.2013.10.008, 2013.

Moecher, D. P. and Samson, S. D.: Differential zircon fertility of source terranes and natural bias in the detrital zircon record: Implications for sedimentary provenance analysis, Earth Planet. Sc. Lett., 247, 252–266, https://doi.org/10.1016/j.epsl.2006.04.035, 2006.

Molnar, P. and England, P.: Late Cenozoic uplift of mountain ranges and global climate change: chicken or egg? Nature, 346, 29–34, https://doi.org/10.1038/346029a0, 1990.

Montgomery, D. R. and Brandon, M. T.: Topographic controls on erosion rates in tectonically active mountain ranges, Earth Planet. Sc. Lett., 201, 481–489, https://doi.org/10.1016/S0012-821X(02)00725-2, 2002.

Oerlemans, J.: Glaciers and climate change, CRC Press, Balkema, Utrecht University, ISBN 902-6-5181-37, 2001.
Östrem, G.: Ice melting under a thin layer of moraine, and the existence of ice cores in moraine ridges, Geogr. Ann., 41, 228–230, https://doi.org/10.1080/20014422.1959.11907953, 1959.

Roe, G. H. and O'Neal, M. A.: The response of glaciers to intrinsic climate variability: observations and models of late-Holocene variations in the Pacific Northwest, J. Glaciol., 55, 839–854, https://doi.org/10.3189/002214309790152438, 2009.

Roering, J. J., Kirchner, J. W., and Dietrich, W. E.: Evidence for nonlinear, diffusive sediment transport on hillslopes and implications for landscape morphology, Water Resour. Res., 35, 853–870, https://doi.org/10.1029/1998WR900090, 1999.

Rowan, A. V., Egholm, D. L., Quincey, D. J., and Glasser, N. F.: Modelling the feedbacks between mass balance, ice flow and debris transport to predict the response to climate change of debris-covered glaciers in the Himalaya, Earth Planet. Sc. Lett., 430, 427–438, https://doi.org/10.1016/j.epsl.2015.09.004, 2015.

Ruhl, K. W. and Hodges, K. V.: The use of detrital mineral cooling ages to evaluate steady state assumptions in active orogens: An example from the central Nepalese Himalaya, Tectonics, 24, TC4015, https://doi.org/10.1029/2004TC001712, 2005.

Schildgen, T. F., Van der Beek, P. A., Sinclair, H. D., and Thiede, R. C.: Spatial correlation bias in late-Cenozoic erosion histories derived from thermochronology, Nature, 559, 89–93, https://doi.org/10.1038/s41586-018-0260-6, 2018.

Schoof, C.: The effect of cavitation on glacier sliding, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 461, 609–627, https://doi.org/10.1098/rspa.2004.1350, 2005.

Schoof, C.: Ice-sheet acceleration driven by melt supply variability, Nature, 468, 803–806, https://doi.org/10.1038/nature09618, 2010.

Sharp, M.: Annual moraine ridges at Skálafellsjökull, south-east Iceland, J. Glaciol., 30, 82–93, https://doi.org/10.3189/S0022143000008522, 1984.

Small, R. J.: Englacial and supraglacial sediment: transport and deposition, in: Glacio-Fluvial Sediment Transfer, edited by: Gurnell, A. M. and Clark, M. J., An Alpine Perspective, Wiley, Chichester, 111–145, 1987.

Spedding, N.: Hydrological controls on sediment transport pathways: implications for debris-covered glaciers, IAHS publication, 133–142, ISBN 190-1-5023-17, 2000.

Steer, P., Huismans, R. S., Valla, P. G., Gac, S., and Herman, F.: Bimodal Plio – Quaternary glacial erosion of fjords and low-relief surfaces in Scandinavia, Nat. Geosci., 5, 635–639, https://doi.org/10.1038/ngeo1549, 2012.

Stock, G. M., Ehlers, T. A., and Farley, K. A.: Where does sediment come from? Quantifying catchment erosion with detrital apatite (U–Th)/He thermochronometry, Geology, 34, 725–728, https://doi.org/10.1130/G22592.1, 2006.

Swift, D. A., Nienow, P. W., and Hoey, T. B.: Basal sediment evacuation by subglacial meltwater: suspended sediment transport from Haut Glacier d'Arolla, Switzerland, Earth Surf. Proc. Land., 30, 867–883, https://doi.org/10.1002/esp.1197, 2005.

Tennant, C., Menounos, B., Ainslie, B., Shea, J., and Jackson, P.: Comparison of modeled and geodetically-derived glacier mass balance for Tiedemann and Klinaklini glaciers, southern Coast Mountains, British Columbia, Canada, Global Planet. Change, 82, 74–85, https://doi.org/10.1016/j.gloplacha.2011.11.004, 2012.

Thomson, S. N., Reiners, P. W., Hemming, S. R., and Gehrels, G. E.: The contribution of glacial erosion to shaping the hidden landscape of East Antarctica, Nat. Geosci., 6, 203–207, https://doi.org/10.1038/ngeo1722, 2013.

Tranel, L. M., Spotila, J. A., Kowalewski, M. J., and Waller, C. M.: Spatial variation of erosion in a small, glaciated basin in the Teton Range, Wyoming, based on detrital apatite (U–Th)/He thermochronology, Basin Res., 23, 571–590, https://doi.org/10.1111/j.1365-2117.2011.00502.x, 2011.

Ugelvig, S. V., Egholm, D. L., Anderson, R. S., and Iverson, N. R.: Glacial Erosion Driven by Variations in Meltwater Drainage, J. Geophys. Res.-Earth, 123, 2863–2877, https://doi.org/10.1029/2018JF004680, 2018.

Ugelvig, S. V., Egholm, D. L., and Iverson, N. R.: Glacial landscape evolution by subglacial quarrying: A multiscale computational approach, J. Geophys. Res.-Earth, 121, 2042–2068, https://doi.org/10.1002/2016JF003960, 2016.

Valla, P. G., Van der Beek, P. A., and Braun, J.: Rethinking low-temperature thermochronology data sampling strategies for quantification of denudation and relief histories: a case study in the French western Alps, Earth Planet. Sc. Lett., 307, 309–322, https://doi.org/10.1016/j.epsl.2011.05.003, 2011a.

Valla, P. G., Shuster, D. L., and Van Der Beek, P. A.: Significant increase in relief of the European Alps during mid-Pleistocene glaciations, Nat. Geosci., 4, 688–692, https://doi.org/10.1038/ngeo1242, 2011b.

Vermeesch, P.: How many grains are needed for a provenance study? Earth Planet. Sc. Lett., 224, 441–451, https://doi.org/10.1016/j.epsl.2004.05.037, 2004.

Whipple, K. X.: The influence of climate on the tectonic evolution of mountain belts, Nat. Geosci., 2, 97–104, https://doi.org/10.1038/ngeo413, 2009.

Willenbring, J. K. and Jerolmack, D. J.: The null hypothesis: globally steady rates of erosion, weathering fluxes and shelf sediment accumulation during Late Cenozoic mountain uplift and glaciation, Terra Nova, 28, 11–18, https://doi.org/10.1111/ter.12185, 2016.

Winkler, S. and Matthews, J. A.: Observations on terminal moraine-ridge formation during recent advances of southern Norwegian glaciers, Geomorphology, 116, 87–106, https://doi.org/10.1016/j geomorph.2009.10.011, 2010.

Yanites, B. J. and Ehlers, T. A.: Intermittent glacial sliding velocities explain variations in long-timescale denudation, Earth Planet. Sc. Lett., 450, 52–61, https://doi.org/10.1016/j.epsl.2016.06.022, 2016.

Zachos, J., Pagani, M., Sloan, L., Thomas, E., and Billups, K.: Trends, rhythms, and aberrations in global climate 65 Ma to present, Science, 292, 686–693, https://doi.org/10.1126/science.1059412, 2001.