Topographic signatures of progressive glacial landscape transformation

The Pleistocene glaciations left a distinct topographic footprint in mountain ranges worldwide. The geometric signature of glacial topography has been quantified in various ways, but the temporal development of landscape metrics has not been traced in a landscape evolution model so far. However, such information is needed to interpret the degree of glacial imprint in terms of the integrated signal of temporal and spatial variations in erosion as a function of glacial occupation time.


Summary
The Pleistocene glaciations left a distinct topographic footprint in mountain ranges worldwide. The geometric signature of glacial topography has been quantified in various ways, but the temporal development of landscape metrics has not been traced in a landscape evolution model so far. However, such information is needed to interpret the degree of glacial imprint in terms of the integrated signal of temporal and spatial variations in erosion as a function of glacial occupation time.
We apply a surface process model for cold-climate conditions to an initially fluvial mountain range. By exploring evolving topographic patterns in model time series, we determine locations where topographic changes reach a maximum and where the initial landscape persists. The signal of glacial erosion, expressed by the overdeepening of valleys and the steepening of valley flanks, develops first at the glacier front and migrates upstream with ongoing glacial erosion. This leads to an increase of mean channel slope and its variance. Above steep flanks and head-walls, however, the observed mean channel slope remains similar to the mean channel slope of the initial fluvial topography. This leads to a characteristic turning point in the channel slopeelevation distribution above the equilibrium line altitude, where a transition from increasing to decreasing channel slope with elevation occurs. We identify this turning point and a high channel slope variance as diagnostic features to quantify glacial imprint.
Such features are abundant in glacially imprinted mid-latitude mountain ranges such as the Eastern Alps. By analysing differently glaciated parts of the mountain range, we observe a decreasing clarity of this diagnostic morphometric property with decreasing glacial occupation. However, catchments of the unglaciated eastern fringe of the Alps also feature turning points in their channel slope-elevation distributions, but in contrast to the glaciated domain, the variance of channel slope is small at all elevation levels.
However, a transient mountain landscape represents just a snapshot of its long-term evolution. The geometry of the original fluvial landscape is in general unknown, hence quantifying the impact of late Cenozoic climate change on mountain topography based on landscape metrics represents a difficult problem. This prevents the interpretation of the evolutionary state of glacial topography. Such knowledge about fluvial landscapes, however, has significantly advanced the field of tectonic geomorphology, and a glacial counterpart is needed. Only a few attempts exist to provide theory (Headley et al., 2012;Prasicek et al., 2020) and applications (Pedersen et al., 2010). Here, we contribute to closing this research gap by reporting the evolution of typical patterns of glacial imprint in the channel slope-elevation domain.
From the perspective of an entire mountain range or even on a global scale, the efficiency of erosional surface processes in cold-versus warm-climate conditions is still debated. While several studies have highlighted a distinct increase in mountain erosion as a consequence of Quaternary climate modulations Molnar, 2004;Zhang et al., 2001), largely uniform erosion rates throughout the Cenozoic are reported by others (Ganti et al., 2016;Koppes & Montgomery, 2009;Larsen & Montgomery, 2012;Willenbring & Jerolmack, 2016). In any case, the transformation from a fluvial to glacial mountain landscape requires strong, valley-scale variations in erosion rate (Champagnac et al., 2014;Montgomery, 2002;Pedersen & Egholm, 2013;Shuster et al., 2011;Steer et al., 2012;Valla et al., 2011). Domains undergoing a fast fluvial-to-glacial transformation can occur in the direct vicinity of persisting fluvial topography. Such 'hot spots' and 'cold spots' in mountain erosion create a characteristic elevation-dependent signature of climate in mountain topography, which results in bimodal landscapes. Landscape bimodality is caused by different patterns of cold-climate erosion above and below the ELA, which will be conceptualized in the following with an emphasis on landscape geometry as captured by the channel slope-elevation distribution (Figure1).
In a fluvial mountain range with increasing uplift rate from the periphery to the centre (tent-shaped uplift) and at (or close to) topographic steady state, the channel slope (i.e. gradient following the steepest descend along the flow paths) and topographic gradient increase monotonically with elevation (Hergarten et al., 2010;Robl et al., 2015aRobl et al., , 2017 (Figure 1 channel slope-elevation plots, dashed line). When climate cools, the ELA drops and cold-climate erosional surface processes progressively transform the fluvial landscape. Findings supporting both the destruction of peak relief and the formation of valley relief allow the definition of two end-member scenarios for cold-climate landscape evolution (e.g. Egholm et al., 2017;Nielsen et al., 2009;Mitchell & Montgomery, 2006;Shuster et al., 2011).
Above the ELA, a combination of glacial and periglacial processes drives the formation of cirques and the destruction of peak relief, resulting in an elevated low-relief landscape towered by steep ridges and horns (e.g. Brozovic et al., 1997;Egholm et al., 2009;Mitchell & Montgomery, 2006;Scherler, 2014) (Figure 1a, the 'glacial buzz-saw').
Hypsometric maxima occur at and above the ELA. The topographic gradient and channel slope follow the original fluvial pattern and increase with elevation up to the ELA, but decrease with the occurrence of cirques above, then increase again in the steep summit domain. This causes a distinct turning point in the channel slope-elevation distribution at, or slightly above, the ELA.
Below the ELA, glaciers transform fluvial valleys into wide, flat glacial troughs. As a result, local base levels are lowered and valley relief is amplified (Figure 1b). This leads to a decrease in topographic gradients (and the channel slope) at valley bottoms, and to a strong increase at the steep valley flanks. A maximum in channel slope occurs due to the strong increase in the channel slope at low-elevation valley flanks relative to the initial landscape. This leads to a turning point in the channel slope-elevation distribution, slightly below the ELA.
In combination, the ratio between the destruction of peak relief and the formation of valley relief controls whether total mountain relief decreases or increases during cold-climate periods (e.g. Robl et al., 2020). With increasing duration of glacial occupation, the combination of both processes ( Figure 1c) should eventually lead to large low-relief areas above the ELA and deep valleys with steep flanks below (e.g. Egholm et al., 2017;Steer et al., 2012). This results in a topographic signal with low mean channel slopes at high elevations but high mean channel slopes at lower elevations (Robl et al., 2015a).
In the following, we systematically explore the progressive transformation from fluvial to glacial landscapes by employing a landscape evolution code capable of simulating glacial and periglacial processes (Egholm et al., 2011. Starting with a fluvial landscape in steady state, we follow the trajectory of landscape change and determine loci and points in time where topographic changes reach a maximum (hot spots) and where the initial landscape persists (cold spots). We trace the fluvial-to-glacial transformation of an entire mountain range and a single valley, and compare derived landscape metrics from model time series with real-world examples of mid-latitude mountain ranges such as the Eastern Alps.

| METHODS
We employ a time-dependent numerical model to describe the evolution of mountain topography during periods of warm-and coldclimate conditions. By applying a set of numerical tools, we analyse modelled and natural mountain landscapes to quantify the impact of cold-climate processes on topography.

| Distribution of channel slope versus surface elevation
We analyse the distribution of channel slope versus surface elevation, an approach that has been applied to both fluvial and glacial terrain (Hergarten et al., 2010;Kühni & Pfiffner, 2001;Robl et al., 2015a). Channel slopes in m/m are calculated by a standard single flow algorithm (D8) implemented in Topotoolbox (Schwanghart & Kuhn, 2010;Schwanghart & Scherler, 2014). All cells defined by a flow accumulation greater than one pixel in the synthetic example and four pixels in the natural landscapes (>0.04 km 2 ) are considered for computing channel slopes. We classify the channel slope data into 50-m-wide elevation bins and report statistical indicators such as mean and percentiles (P25, P50 and P75) for each bin. For a graphical representation of the channel slope distribution per elevation bin, the relative frequency of the occurrence of channel slope values is evaluated using channel slope bins with a width of 0.02 m/m (Robl et al., 2015a). Further, we investigate the shape of the entire channel slope-elevation distribution.

| Landscape evolution model
The computational landscape evolution model employed here can simulate glacial, fluvial and hillslope erosion. The fluvial component is only used here to provide a fluvial steady-state topography as a reference and starting point to systematically explore the fluvial-to-glacial landscape transformation. To obtain a fluvial steady-state landscape, we employ a detachment-limited model for bedrock channel incision (e.g. Braun & Willett, 2013;Hergarten, 2002;Howard, 1994Howard, , 1980. Erosion in the hillslope domain is approximated by diffusion to produce the initial fluvial steady state. In the glacial model, the flow of ice is handled by an integrated second-order shallow ice approximation scheme (iSOSIA) (Egholm et al., 2011). This approach aims to balance the requirements of accuracy in simulating the ice flow and computational efficiency in modelling the long-term landscape evolution. In contrast to standard shallow ice approximation (SIA) models, iSOSIA considers longitudinal and transverse stresses in the ice and hence resolves horizontal stress gradients, which are important for alpine settings, because they provide local coupling of ice flow and impede strong local gradients in ice velocity (Egholm et al., 2012b). Accumulation and ablation of ice are modelled by a simple positive-degree-day approach, and additionally modified by avalanches and snow drift . We model ice temperate and ignore the effects of basal freezing at the upper reaches of thin, high-altitude alpine glaciers. Glacial erosion processes are represented by abrasion and quarrying. We model glacial erosion on bedrock and neglect the potential effects of subglacial sediment protecting the glacier bed. Abrasion is estimated based on ice sliding at the glacier base (Egholm et al., 2012a;Hallet, 1979;Lliboutry, 1994). We use a sliding law that accounts for the opening of cavities due to the roughness of the bed (Ugelvig et al., 2018b).
Consequently, subglacial hydrology controls the effective pressure of F I G U R E 1 Sketch of the two endmember models (a, b) and the combined model (c) on the evolution of bimodal landscapes explained in the text. Red circles indicate the turning point from increasing to decreasing channel slope, which mark the transition from less to high effective glacial erosion (a) and vice versa (b) [Color figure can be viewed at wileyonlinelibrary.com] the system and influences the basal sliding speed . Subglacial quarrying is additionally influenced by effective pressure directly and bed channel slope in the direction of sliding (Iverson, 2012;Ugelvig et al., 2016). In the glacial model, a temperature-dependent frost cracking approach is employed to model periglacial erosion Egholm et al., 2015). Since steep surface slopes dramatically reduce the performance of the numerical model, a local slope threshold criterion is used. This criterion immediately lets all local slopes steeper than a threshold of 1 (45 ) collapse to the threshold slope.
While this criterion is required for numerical reasons, it can be interpreted as some kind of landsliding.
Tectonic uplift is only considered to produce the initial fluvial topography, while isostatic effects related to loading/unloading of the lithosphere due to ice and erosion are also included.
Detailed description of the combined model can be found elsewhere (Braedstrup et al., 2016;Egholm et al., 2011Egholm et al., , 2012aEgholm et al., , 2015Ugelvig & Egholm, 2018a). However, a summary of the setup and chosen parameter values can be found in the Supplementary Information. We use standard parameter sets for all model components. In the following, we briefly describe the most important aspects of the model setup.

| GENERAL MODEL SETUP
The initial topography for all the synthetic experiments has an extent of 60 Â 80 km and consists of 10,000 irregular Voronoi cells. We start with a 500-m elevated plain and some low-amplitude random noise. we set a constant sea level temperature of 6.5 C and 8.5 C, respectively, a lapse rate of 6 C/km, a seasonal temperature variation of 14 C and a positive-degree-day melt rate constant of (0.001 m d À1 C À1 ) , which results in an ELA of $1300 m and $1650 m. Precipitation is set to 1 m/a for the entire model domain. These mass-balance parameters are chosen to cover the landscape with an extended ice layer so that the initial ice accumulation areas in both experiments has roughly the same size. The free-scaling erosion parameters K q and K a are uniform in all our experiments as we are interested in glacial erosion patterns controlled by pre-glacial topography and not on the impact of substrate properties (i.e. lithology). They are set to 1 Â 10 À3 Pa À1 and 1 Â 10 À4 Pa À1 , respectively.
Our model setup provides little capacity for frost cracking, because a constant mean annual temperature limits the landscape to be exposed to frequent freeze-thaw cycles (Anderson, 1998;Hales & Roering, 2007). In addition, a thick ice cover dampens the surface temperature variations needed for promoting segregation ice growth.
To account for the glacial-interglacial cycles with waxing and waning glaciers typically observed in glaciated mountain ranges at mid-latitudes, we perform a third experiment (c) where we introduce three asymmetric, sawtooth-shaped temperature cycles with a period of 100 ka and an amplitude of 3 C. Each cycle is represented by 80 ka of constant cooling followed by 20 ka of rapid warming, which results in an ELA varying between $1300 and 1850 m. We stress that this approach is not meant to simulate a specific time period but simply to account for general temperature variations in the long-term climate record (i.e. Lisiecki & Raymo, 2005) and their specific effects on topography (e.g. cirque formation).

| RESULTS
In the following, we explore the emerging spatial and temporal patterns and changes in landscape geometry. First, we highlight topographic changes due to glacial and periglacial erosion on a mountain range scale. Second, we track temporal and spatial variations in glacial erosion and the resulting topographic patterns in a single catchment.
Finally, we compare the model results with selected catchments of the Eastern Alps to estimate the degree of glacial reshaping.

| Spatio-temporal pattern of glacial erosion on a mountain range scale
The initial fluvial mountain range is dominated by a west-easttrending principal drainage divide with a maximum elevation of  Once cold-climate surface processes start reshaping fluvial topography, the mean and variance of the channel slope increase above the ELA (Figure 4b). This trend is amplified with ongoing glacial occupation. After 300 ka of glacial erosion, the topography is represented by a channel slope-elevation distribution with a distinct channel slope maximum slightly above the ELA and a second maximum below the initial base level where glacial troughs formed (Figure 4c, d).
The spread between the first and third channel slope quartile distinctly increases from the initial fluvial to final glacial landscape across all elevations, with the strongest rise at and above the ELA (Figure 4c).
In addition to the general channel slope-elevation distribution, In general, the pattern of glacial erosion shows strong spatial variations. This causes a significant increase in local relief from the original fluvial to the resulting glacial landscape (Figure 5e and f, respectively).

| The signature of glacial erosion in three alpine catchments
We compare the progressively evolving glacial signature from our modelled synthetic landscapes with landscape characteristics from catchments of the Eastern Alps to determine the degree of glacial imprint in mid-latitude mountain ranges. In the Eastern Alps, fully, partly and non-glaciated catchments with similar peak elevation, substrate properties and structural inventory exist in spatial proximity ( Figure 8a). This allows elimination of the otherwise strong lithological control on mountain topography (i.e. the persistence of fluvial and glacial landforms) (e.g. Prasicek et al., 2015;Robl et al., 2015a).
We choose the Trofaiach catchment (Mur) (Figure 8c Anderson et al., 2006;Braedstrup et al., 2016;Egholm et al., 2011;Harbor et al., 1988). Here, we discuss the spatial and temporal patterns of glacial erosion emerging from our numerical experiments, the characteristics of the resulting synthetic landscapes and channel slope-elevation plots to diagnose the evolutionary state of real-world topography.

| Spatial variations in erosion rate
The spatial pattern of glacial erosion is controlled by climate, i.e. the ELA and the mass-balance gradient with elevation, and the initial fluvial landscape. However, the combined influence of these factors can cause an increase (Montgomery, 2002;Shuster et al., 2011;Valla et al., 2011) or decrease (Brozovic et al., 1997;Mitchell & Montgomery, 2006) of relief. Both effects have been observed in high-latitude landscapes, with an increase of fjord relief and a decrease of interfluve relief above (Egholm et al., 2017). Our simulations inspired by mid-latitude mountain ranges dominated by alpine-style glaciations with peaks and ridges towering above a warmbased glacier network show a general increase in relief.
We observe selective glacial erosion where fluvially predefined channels are excavated faster and deeper than tributaries and the surrounding topography (Figures 6 and 3). In our simulations, this results from the spatial distribution of ice flux, basal water availability and basal glacier velocity. Domains along the main channels undergo a strong fluvial-to-glacial transformation, while areas directly above with low ice flux maintain their fluvial topographic signature and hanging sections form at tributary intersections.
The vertical incision in valleys causes lengthening of flanks and headwalls and eventually forms persisting peaks that are surrounded by steep surfaces, which promote snow avalanching and reduce glacial erosion (Ward et al., 2012). However, if these ice-free flanks are within the elevation range to enable frost cracking, more rapid headwall retreat by undercutting should be possible, as the glacier below would transport the debris away and prevent the formation of scree deposits (Scherler, 2014). In general, the ratio between glacial and periglacial erosion may differ depending on the conditions at the glacier base (Egholm et al., 2017) and efficiency of frost cracking in ice-free areas (Hales & Roering, 2007;Scherler, 2014). In our simulations, the efficiency of periglacial erosion is higher under varying than constant climate conditions. During interglacial periods, the glacially steepened valley sides experience enhanced frost cracking after deglaciation. However, in both scenarios, periglacial erosion amounts to only a few metres and hence is negligible compared with the amount of glacial incision amplifying relief.
The widening and deepening of valleys in concert with the preservation of summit domains decrease the mean elevation and reduce the topographic load of mountain ranges. Isostatic compensation of localized erosion in valleys can enhance this effect and raise surrounding peaks .

| Temporal variations in erosion rate
The rates of topographic adjustment to glacial conditions vary in space and time. With ongoing glacial occupation, erosional hotspots with erosion rates distinctly exceeding background rates evolve and progressively form overdeepenings, which grow both in depth and laterally.
Glacial erosion rates have been constrained in general (e.g. Cook et al., 2020), but temporal variations in erosion rate showing the development and migration of erosional hotspots have only been documented in rare cases of extraordinary deep glacial valley incision (Haeuselmann et al., 2007;Shuster et al., 2011;Valla et al., 2011). For most parts of glacial landscapes, temporal changes in erosion rate remain hardly constrained.
Our numerical experiments contribute to solving this issue.
As the complexity of glacial erosion processes is poorly understood, glacial landscape evolution models generally use simple erosion rules that relate the efficiency of erosion to sliding at the glacier base (Hallet, 1979;Harbor, 1992;Herman et al., 2011). Besides all their simplifications and limitations, such models allow constraining the timing and rate of glacial landscape adjustment (e.g. Pedersen & Egholm, 2013;Sternai et al., 2013). In our experiments, quarrying and abrasion represent the primary mechanisms for subglacial erosion, whereby basal sliding speed and hence the related basal shear stress are primary factors of their efficiency and control the spatio-temporal pattern of topographic adjustment from a fluvial to glacial landscape.
Our experiments along with former numerical modelling studies demonstrate that subglacial shear stress decreases as glacial erosion progressively transforms preglacial V-shaped valleys into U-shaped troughs (Braedstrup et al., 2016). This is also consistent with our findings showing that the rate of channel slope adjustment is distinctly higher in the first 100 ka of glacial occupation than in the last 200 ka ( Figure 4e and Supporting Information Movie S2). As a result, mean erosion rates are highest in the early stages of glacial landscape development, when the topography is not yet adapted to the new cold-climate conditions (Harbor, 1992). Spatial variations in erosion rate and the timing of maximum rates are dominated by valley overdeepening (Figures 3,6 and 7). Overdeepening first occurs in down-glacier reaches and progressively extends upstream (Figure 7b).
We observe a strong positive feedback between the formation of overdeepenings and the sliding rate, indicating that valley excavation by glaciers is a non-linear process. Once overdeepening is initiated, the valley floor is rapidly lowered (Figure 7b, blue, magenta and red curve). Glacial hydrology plays a key role in this behaviour, as it affects subglacial sliding and thereby influences long-term landscape evolution. The highest volumes of meltwater are available in the lower part of the glacial system, increasing sliding rates and leading to a boost in erosion . Lowering of the glacier bed also lowers the glacier surface and hence leads to increased melting, again increasing sliding. Starting at the terminus, this signal of glacial erosion migrates upwards while the glacier progressively eats away its own base ( Figure 7b) and retreats. This mechanism has the double effect of decreasing erosion rates at low elevations and triggering headward propagation of glacial erosion, which has also been shown by former glacial modelling (Sternai et al., 2013) and erosion studies . Therefore, our results suggest that glacial valley deepening and erosion at and above the ELA are related processes.
Assuming constant cold-climate conditions, our model predicts only Further, headwall retreat occurs to a large part above the ice, where periglacial processes dominate erosion. This was predicted by numerical modelling (e.g. Anderson et al., 2013;Egholm et al., 2015;Hales & Roering, 2007), and shown by cosmogenic nuclide (Mair et al., 2020) and monitoring studies (Hartmeyer et al., 2020). Thus, the few metres of periglacial erosion predicted in our experiments might strongly underestimate real-world headwall retreat, at least in temperate mountain ranges. Further, the hypsometry of the initial landscape controls the mass balance and the spatial distribution of the glacial erosion potential. In a steady-state fluvial mountain range, this potential is small at highest elevations consisting of steep ridges and peaks, while it can be large in a premature landscape with elevated plateaus (e.g. Pedersen et al., 2014).

| Channel slope-elevation distributions for interpreting the degree of glacial imprint
In our experiments, the glacial landscape geometry is finally characterized by spacious, flat valleys separated from persisting summit areas by (high) steep flanks (Figures 6 and 3). The development of such neighbouring erosional 'hot spots' and 'cold spots' creates a characteristic signature in the evolving glacially conditioned mountain range. Such a signature is visible in the channel slope-elevation distributions of the synthetic and the investigated natural mountain landscapes. In particular, it is expressed by (Figure 4e Thus, the inflection point in the channel slope-elevation distribution around the ELA is rather due to reducing the proportional fraction of topography below and above than a net accumulation of flat surfaces around. Mountain ranges worldwide (e.g. European Alps, Andes, Tien Shan, Pyrenees) are characterized by a similar transition from increasing to decreasing channel slopes with elevation, while channel slopes increase continuously in others (e.g. Atlas, Carpathians) (for review Robl et al., 2017). In addition to the glacial steepening of low-lying landscape parts modelled and observed in our study (Figures 4 and 9), other influences have been suggested to cause such contrasting channel slope distributions (e.g. Hergarten et al., 2010;Kühni & Pfiffner, 2001;Robl et al., 2015a).
According to the hypothesis of fluvial prematurity, low-gradient high-altitude surfaces are remnants of low-relief topography lifted by a recent uplift event (e.g. Hack, 1975;Hergarten et al., 2010). In this case, gentle channel slopes and valleys form the summit domains, while the regions below are steep and incised, as the two contrasting landscape types represent the ancient and recent tectonic regime, respectively. This type of elevated low-relief surfaces is not related to the glacial extent, shows a lower interquartile range than glacial landscapes in the channel slope-elevation distribution and lacks steep summit domains.
The glacial buzzsaw hypothesis on the other hand assumes longterm glacial erosion to decrease relief and channel slope above the equilibrium line altitude, which eventually results in a concentration of surface area at this elevation (Figure 1b) (e.g. Brozovic et al., 1997;Egholm et al., 2017;Mitchell & Montgomery, 2006). This type of elevated low-relief landscape occurs within a glaciated part of the mountain range, relates to the ELA (at and above) and shows a very high interquartile channel slope range. The channel slope-elevation distributions thus serve to distinguish between different states and processes dominating landscape evolution.

| Channel slope-elevation distributions of the Eastern Alps
The channel slope-elevation characteristics should in principle allow the interpretation of the degree of glacial imprint of entire mountain ranges. During the LGM, the Eastern Alps featured a continuous transition from a fully glaciated western to a fully unglaciated eastern part.
The varying degree of glacial imprint is deducible today by analysing the channel slope-elevation distributions of individual catchments.
The three catchments analysed in this study, representing the minorly However, a similar pattern but with a significantly lower interquartile channel slope range also emerges at the eastern fringe of the Alps outside the LGM ice extent (Robl et al., 2015a). This advocates the interpretation of a premature landscape induced by a tectonic driver (Hergarten et al., 2010;Legrain et al., 2014;Robl et al., 2015a;Wagner et al., 2010). In this part of the Eastern Alps, the vertical position of the turning point in the channel slope-elevation distribution varies. In contrast, the glaciated realm of the Eastern Alps shows significantly higher interquartile channel slope ranges and the turning point, generally located at about 1800 m (Figure 8h), roughly marks the altitude of the ELA (Robl et al., 2015a). This highlights the diagnostic capabilities of channel slope-elevation distributions.

| Limitations
To study the feedbacks between glacial dynamics and topographical development, we kept our experiments as simple as possible and neglect several processes that influence landscape evolution. For example, fluvial erosion where glaciers are absent (e.g. Schlunegger & Norton, 2013), tectonic forcing (e.g. Sternai et al., 2013;Tomkin et al., 2003), stabilizing feedbacks in modelling erosion for adverse slopes within overdeepenings (e.g. Alley et al., 2003;Haeberli et al., 2016;Patton et al., 2016;Werder, 2016) and subglacial bedrock fracturing in response to high differential stresses (e.g. Leith et al., 2014) are not included.
Further, the erosion laws' empirically determined proportionality constants and velocity exponents are particularly uncertain, ranging roughly from 0.6 (Cook et al., 2020) to 2 (e.g. Herman et al., 2015;Koppes et al., 2015). However, variations in parameters such as erosion law scaling constants only control the pace of erosion, while the spatial patterns of erosion are qualitatively similar (Braedstrup et al., 2016;Egholm et al., 2011Egholm et al., , 2012a. Therefore, we focus our analysis on the spatial and temporal patterns rather than absolute rates of erosion. Hillslope diffusion is only included to maintain maximum bed channel slopes of 45 , which is considered the upper limit for iSOSIA's validity (Egholm et al., 2011). This has a direct effect on the landscapes' geometry. With continuously ongoing glacial erosion on the valley floor and persisting topography above, surface gradients are rising and the valley flanks tend to steepen beyond the slope threshold. As a consequence, hillslope processes are turned on to prevent the landscape from becoming too steep, and hence widen the valley with uniform flank side gradients of 45 . This results in channel slopes clustering around 45 within all altitude ranges ( Figure 4) and a slight conversion of emerging U-shaped valleys towards more V-shaped valley crosssections (Figures 3 and 6). However, the final modelled landscape (Figure 2b) possesses only a small surface fraction with bed channel slopes of 45 , and allowing the topography to rise above this threshold channel slope would not condemn the diagnostic capabilities of channel slope-elevation distributions and may even clarify the distinction between highly and less glacially modified landscape patches.

| CONCLUSIONS
Numerical modelling of glacial and periglacial erosion processes on a synthetic and a natural fluvially conditioned landscape combined with morphological analysis of topographic properties of the Eastern Alps leads us to the following conclusions: The model results show that glacial erosion is focused in valleys of the fluvially conditioned drainage system, with the highest total erosion at trunk valleys (erosion hot spots). Peaks and ridges towering above the glacier network are only affected by periglacial erosion, which amounts to only a few metres (erosion cold spots). Thus, we observe a general increase in relief due to glaciation in these regions.
Spatial variations in erosion rate and the timing of maximum rates are dominated by valley overdeepening. Overdeepening first occurs in down-glacier reaches and progressively extends upstream. Controlled by the hydrological configuration of the glacier base, we observe a strong positive feedback between the formation of overdeepenings and the sliding rate. This makes valley excavation by glaciers a nonlinear process.
In our simulation, glacial erosion leads to a highly bimodal landscape with the highest channel slopes at steep valley flanks and tributary headwalls, whereas areas above maintain their fluvial characteristics with rather low channel slopes gradients. Two transitions from increasing to decreasing channel slopes with elevation are emerging, one below (where overdeepenings dominate the channel slope distribution) and one above the ELA. In natural landscapes, the first turning point is rarely observable as overdeepenings are filled by water and sediment.
Glacial processes alter the (large-scale) topographic pattern of initial fluvially conditioned landscapes. Distributions of channel slope with elevation vary with the degree of glacial imprint, as demonstrated by the progressive evolution from a fluvial to glacial landscape in our numerical experiments. This temporal transition can also be observed in spatial vicinity from a fully glaciated western to a minorly glaciated eastern part of the Eastern Alps. Thus, channel slope-elevation plots serve as a diagnostic tool for interpreting the glacial-fluvial influence in mountain landscapes.
The position of the turning point and the variance of channel slopes in the channel slope-elevation plots serve to distinguish between drivers of the formation of bimodal landscapes. In a glacial landscape, the turning point relates to the ELA and shows a high interquartile channel slope range. On the contrary, in bimodal landscapes without glacial imprint, the vertical position of the turning point varies and the interquartile range is comparatively small.
While the topographic evolution of fluvially conditioned mountain ranges has been successfully explored by the application of stream power incision models (SPIMs) (e.g. Robl et al., 2017, and references therein), a high-resolution, numerical description of glacial mountain ranges and their evolution towards a hypothetical glacial steady state is still in its infancy. In contrast to the SPIMs, where fluvial erosion is computed from geometrical properties of the drainage system (i.e. local channel slope and upstream contributing drainage area), state-of-the-art glacial erosion models are computationally expensive as they describe erosion based on the physical process of moving ice over its bedrock (e.g. Braun et al., 1999;Egholm et al., 2011;Herman et al., 2011). These models face some limitations and stability issues (e.g. steep slopes) when considering the mathematical and numerical complexity of ice dynamics. While these models are extremely useful for describing landscape modifications in periods of fluctuating climate (e.g. during a glacial-interglacial cycle), computational performance and numerical stability impede high-resolution, large-scale simulations over the million years' time scale. Recently and similar to the SPIM, a new approach to understanding glacial landscape evolution has been published (Deal & Prasicek, 2021), which allows the computation of glacial erosion rates from topographic properties and simple climate assumptions. Simultaneously, Hergarten (2021) presented a new landscape evolution model, describing glacial landscape evolution based on a stream power law for glacial erosion, and coupled this model with a fluvial erosion model taking into account both incision and sediment transport. Further benchmarking will show whether these models are able to produce large-scale landscape patterns observed in natural glacially shaped landscapes and readily described by higher-order ice flow models for small domains. However, if these rather simple models are similarly successful as their fluvial counterparts, we will apply this type of models and build on the findings of this study to explore the processes and topographic signatures during the transition from fluvial to glacial landscapes (and vice versa) for large mountain ranges with high spatial resolution over the million years' time scale.

DATA AVAILABILITY STATEMENT
The code (iSOSIA) and the experimental setup is available on request from D.L.E. and M.L.