The periglacial engine of mountain erosion – Part 2 : Modelling large-scale landscape evolution

There is growing recognition of strong periglacial control on bedrock erosion in mountain landscapes, including the shaping of low-relief surfaces at high elevations (summit flats). But, as yet, the hypothesis that frost action was crucial to the assumed Late Cenozoic rise in erosion rates remains compelling and untested. Here we present a landscape evolution model incorporating two key periglacial processes – regolith production via frost cracking and sediment transport via frost creep – which together are harnessed to variations in temperature and the evolving thickness of sediment cover. Our computational experiments time-integrate the contribution of frost action to shaping mountain topography over million-year timescales, with the primary and highly reproducible outcome being the development of flattish or gently convex summit flats. A simple scaling of temperature to marine δO records spanning the past 14 Myr indicates that the highest summit flats in midto high-latitude mountains may have formed via frost action prior to the Quaternary. We suggest that deep cooling in the Quaternary accelerated mechanical weathering globally by significantly expanding the area subject to frost. Further, the inclusion of subglacial erosion alongside periglacial processes in our computational experiments points to alpine glaciers increasing the long-term efficiency of frost-driven erosion by steepening hillslopes.


Introduction
The question of whether late Cenozoic global cooling caused accelerated erosion of Earth's mountain landscapes is widely debated (Zhang et al., 2001;Molnar, 2004;Willenbring and von Blanckenburg, 2010;Herman et al., 2013;Von Blanckenburg et al., 2015).Support for accelerated erosion comes primarily from ocean basin sedimentary records (Hay et al., 1988;Zhang et al., 2001;Molnar, 2004) and from thermochronological studies of mountainous areas (Herman et al., 2013).Both sources of information indicate the onset of accelerated erosion around 2-6 Myr ago, although resolving long-term transients in both sediment deposition and bedrock erosion is hampered by methodological challenges (e.g.Gardner et al., 1987;Sadler and Jerolmack, 2015), The mechanisms that might drive accelerated erosion are not well understood, but the increase is most clearly detected in glaciated regions, leading many to suggest that global cooling enhanced the capability of glaciers to carve valleys and cirque basins (Haeuselmann et al., 2007;Valla et al., 2011;Herman et al., 2013).The role of glaciations is backed by a close temporal correspondence between the estimated total erosion and the δ 18 O curves that serve as a proxy for palaeotemperature and land-ice volume (Herman et al., 2013).Other authors argue that the main driver of accelerated erosion is the increased amplitude and frequency of climatic oscillations that accompany the cooling trend (e.g.Zhang et al., 2001;Molnar, 2004).A high level of climatic variability can lead to a situation in which landscapes are rapidly readjusting to new boundary conditions and dominant erosional processes (e.g.switching between glacial and fluvial erosion), which may result in persistently high erosion rates (Braun et al., 1999).
The notion that cold-climate processes are responsible, in large part, for accelerated late Cenozoic erosion is supported by the topography of mountain ranges.Most high mountain ranges were glaciated at some time in the Quaternary, and glaciers are thought to have reshaped landscapes significantly within a few million years (e.g.Penck, 1905;Harbor et al., 1988;Hallet et al., 1996).Fjords, U-shaped troughs, hanging valleys, and cirques represent distinct products of efficient subglacial erosion, and these landforms are now prominent in many mid-to high-latitude mountain ranges (e.g. in Norway, Greenland, British Columbia, New Zealand, and southern Chile).Nevertheless, two general observations suggest that processes not directly associated with glaciers could also be of importance for accelerated erosion: (1) a large proportion of low-and mid-latitude mountain ranges have never hosted large ice masses or bore them only briefly and (2) global cooling over the last 15 Myr (Zachos et al., 2001) caused a vast areal expansion of land surface experiencing frost-driven weathering and sediment transport.Under the current warm interglacial conditions, ∼ 18 % of Earth's icefree land surface has a mean annual air temperature below 0 • C (calculated from data in Hijmans et al., 2005), representing a huge potential source of frost-driven sediment production and supply.Consequently, we suggest that the contribution of periglacial processes to global sediment flux over the late Cenozoic may be notably underestimated.
Measurements of millennial-scale erosion rates (Delunel et al., 2010) as well as documented links between elevation and rock scree distribution (Hales and Roering, 2005), weathering zones (Savi et al., 2015), and threshold slopes (Scherler, 2014) suggest that frost-driven erosion can be efficient (on the order of 1 mm a −1 ) in steep mountain landscapes.Frost-related processes may therefore conceivably exert first-order control on the evolution of steep glacial headwalls and valley sides in many areas.
Cold-climate, non-glacial surface processes are also associated with widespread areas of low-relief at high elevations known as summit flats (Fig. 1), which are characteristic of alpine landscapes in, for example, the Laramide ranges of the western USA (e.g.Small and Anderson, 1998;Anderson, 2002;Anderson et al., 2006;Munroe, 2006) and the Caledonian mountains in Scotland, Scandinavia, and Greenland (e.g.Rea et al., 1996;Fabel et al., 2002;Phillips et al., 2006;Nielsen et al., 2009;Goodfellow, 2012).Summit flats are typically mantled thinly with regolith that is widely held to be the product of frost-driven physical weathering (Anderson, 2002;Ballantyne, 2010;Goodfellow, 2012), although there is disagreement concerning the timing of their formation, the contribution of chemical weathering, and relationships to climate (Rea et al., 1996;Whalley et al., 2004;Strømsøe and Paasche, 2011;Goodfellow, 2012).
Previous studies have shown that models of regolith production and transport by frost can explain the regolith cover as well as the smooth, convex shape of the summit flats (Anderson, 2002;Anderson et al., 2012).However, substantial questions remain concerning the quantities of sediment produced and the role of these high summit flats in feeding debris to paraglacial sediment systems in the valleys below (e.g.Ballantyne, 2010).Another aspect concerning geomorphic activity is that high "accordant surfaces" composed of many separate summit flats at a similar elevation are frequently used as indicators of discrete tectonic uplift events (e.g.Widdowson, 1997;Lidmar-Bergström et al., 2000;Bonow et al., 2003Bonow et al., , 2006;;Japsen et al., 2009).The basic assumption behind this "tectonic hypothesis" is that low-relief landforms can form only at base (sea) level, and therefore the presence of low-relief topography at high elevations is best explained by vertical displacement of the landscape after it was formed (Lidmar-Bergström et al., 2013).These unresolved questions further motivate a better understanding of the nature of summit flats because if low-relief landforms can develop in situ at high elevations via frost action, this must be taken into account before uplift is inferred from topography alone.
Analysis of cosmogenic nuclide concentrations in bedrock and boulders from cold-region summit flats has revealed that (1) surface exposure ages often predate considerably the most recent Pleistocene glaciation, and (2) weathering processes must be very slow, with erosion rates not more than a few tens of m Myr −1 (e.g.Small et al., 1997;Fabel et al., 2002;Phillips et al., 2006;Goodfellow et al., 2014).The summit flats are therefore generally considered to be slowly evolving elements within landscapes that juxtapose deeply incised valley troughs cut by Pleistocene glaciers (Small and Anderson, 1998;Anderson, 2002).Understanding the long-term topographic development of these enigmatic landforms remains a challenge.Recent studies of the Norwegian mountains suggest that, although characterised by low erosional activity during recent glacial cycles, high-elevation, low-relief areas may have contributed notable sediment volumes to offshore basins over the last 2.8 Myr (Steer et al., 2012).It seems that accounting for the long-term development of the summit flats requires the processes at play to be integrated throughout the Quaternary period, and possibly even further back because periglacial erosion may have been active several million years before glaciers were introduced into the landscape (Nielsen et al., 2009).
The aim of this study is to time-integrate frost-driven processes and explore landscape-scale feedbacks among frostweathering intensity, sediment mobility, and the evolution of relief in response to changing climate.To achieve this we implement a frost-cracking and frost-creep model into a largescale landscape evolution model.We devise three experimental scenarios that illustrate a range of topographic and climatic boundary conditions reflecting natural landscapes: (1) a small-scale (2 km × 3 km grid), low-relief landscape resembling an individual summit flat surface; (2) a larger-scale (25 km × 50 km grid), high-relief landscape subject to a climate forcing that resembles late Cenozoic cooling; and (3) a large-scale (25 km × 50 km) landscape where periglacial and glacial erosion processes operate in tandem.

The landscape evolution model
The computational landscape evolution model operates on a staggered regular grid of cells (Fig. 2).The bedrock elevation, b(x, y), and the sediment thickness, S(x, y), are recorded in the midpoint of each cell.The surface elevation, h(x, y), of a grid cell is then We use the grid structure to model the evolution of sediment thickness and bedrock elevation over time.The term sediment as used here equates with mobile regolith; that is, unconsolidated minerogenic material able to be mobilised by surface transport processes.
(i, j) The periglacial landscape evolution model applies a staggered grid in which bed elevation, b, sediment thickness, S, and weathering rate, ė, are stored in cell midpoints.The horizontal sediment flux components, q x s and q y s , are computed at cell edges.
The continuity equation quantifies changes in sediment thickness over time.ė is the soil production rate (the rate of bedrock conversion to sediment), ρ s = 2400 kg m −3 is the prescribed (average) bulk density of sediment, and ρ r = 2900 kg m −3 is the density of bedrock.qx and qy represent the sediment flux along the two horizontal coordinate axes x and y.The sediment flux components are calculated at the edges between cells (Fig. 2).We use the frost-cracking intensity (FCI) and the sediment transport efficiency (κ), described in the companion paper (Andersen et al., 2015), to compute the weathering rate and sediment flux at every grid point: ė(x, y) = k e FCI(x, y) (3) We assume that the weathering rate scales linearly with the frost-cracking rate, and k e is the free scaling parameter.q s is the sediment flux parallel to the surface in the down-slope di- Frost-cracking intensity (FCI) Frost transport efficiency (κ) a b where θ = tan −1 ( ∂h ∂x ) and ψ = tan −1 ( ∂h ∂y ).As documented by Andersen et al. (2015), both FCI and κ are estimated for a wide range of mean annual temperatures (MATs) and sediment thicknesses (S) (Fig. 3).Other parameters, such as the thermal properties of sediment and bedrock, the amplitudes of annual and daily temperature variations, as well as the availability of water, also affect FCI and κ, but we do not pursue those issues here.We focus instead on how feedbacks between landscape evolution and frost-driven erosion operate across temporal and spatial variations in temperature and sediment thickness.In every time step of a model simulation, we compute, by interpolation of the rates shown in Fig. 3, the FCI for each cell midpoint and κ for every cell edge using the MAT and sediment thickness, S, of the grid points.
The sediment thickness of each cell is thereafter updated using explicit time integration, and the elevation of the bedrock surface is lowered by Due to numerical stability requirements, the explicit time integration limits the length of the time step to where min( x, y) is the minimum grid spacing and max(κ) is the maximum sediment diffusivity.
In addition to this time-step constraint, the sediment flux is limited by the sediment available in the cells.While conserving the sediment volume, this constraint prevents the occurrence of negative sediment thickness values.As an example, we write the flux constraints for the size of the horizontal flux qx i,j + 1 2 i,j x/ t for ∂h/∂x < 0 S t i,j +1 x/ t for ∂h/∂x > 0 . (10)

Computational experiments
We report the results of three types of computational experiments that all apply the expressions for frost cracking and frost creep (Fig. 3) in order to integrate erosion and sediment transport over time.The experiments are designed to explore how the combined effects of frost cracking and frost creep influence landscape evolution over several million years in the absence of other surface processes; only in the third experiment do we go beyond the periglacial system to include glacial erosion.We acknowledge the contribution of other surface processes to the evolution of landscapes, such as fluvial and aeolian erosion and transport, mass wasting, and chemical weathering; however, we explicitly omit them from our experiments in order to increase transparency.Likewise, we do not consider complicating factors such as differences in rock mass strength that arise via different lithologies and tectonic deformation histories.We strive to apply a firstprinciples approach to study (1) the topographic fingerprint of frost cracking and frost creep in relation to temperature variations that mimic those of the late Cenozoic, (2) the factors limiting long-term erosion rates, and (3) how late Cenozoic temperature shifts may have affected sediment production by frost cracking in realistic landscapes.

Experiment 1: the evolution of periglacial summit flats
The first experiment illustrates how the periglacial frost activity represented by Fig. 3 leads naturally to smooth parabolic areas capped by a thin sediment cover.According to Anderson (2002), this style of periglacial landscape evolution decreases local relief and establishes a steady-state landform, which thereafter is uniformly lowered by slow erosion.Through several model simulations, we first explore how weathering rate and the emergent sediment thickness depend on the climatic conditions and the free parameter k e that scales the bedrock weathering rates.We initiate the first simulation using a fluvial-style landscape in a 2 km×3 km grid (i.e.100×150 cells with a spatial resolution of 20 m; Fig. 4a).The initial relief was generated by a landscape evolution model based on a fluvial streampower erosion law (Braun and Sambridge, 1997).It includes a narrow central ridge and valley structures that reach from the grid centre line to the left and right grid boundaries.The maximum initial relief is 200 m and the mean slope is 20 %, which represents relatively rugged topography.The morphology of the initial relief is not of particular importance, provided that it differs from the smooth parabolic surfaces that represent the model output.We simulate frost cracking and frost creep during recurrent climate cycles that last 100 kyr and entail linear MAT variations between −6 and 0 • C. In order to keep this experiment as simple as possible, temperature is not varied across the surface but through time only.The total simulation period is 4 Myr (i.e.40 climate cycles).During the simulation the sediment thickness is forced to 0 at the left and right boundaries, resembling topographic edges where the sediment from the modelled surface drops into deeper valleys.This type of boundary condition is particularly relevant for high surfaces bounded by steep glacial troughs (Fig. 1a).For the first simulation k e = 10 −3 • C −1 a −1 .This value for k e is not constrained by measurements, and we therefore vary the parameter in subsequent experiments in order to study its influence.
During the simulation, the original topography of the surface is erased, leaving a smooth low-relief convex surface (Fig. 4).Initially, frost cracking attacks the high ridges while the associated debris accumulates in the small valleys between the ridges.Frost cracking is limited in valley floors owing to the thick sediment cover, and topographic relief drops quickly as erosion of the bedrock ridges continues (Fig. 5).After almost 4 Myr (40 climate cycles), the landscape reaches a steady state where sediment thickness (∼ 3 m), relief, and surface curvature are nearly uniform across the landscape (Fig. 5).This configuration is independent of the initial topography and represents the exclusive outcome of the simulated periglacial processes.These results are in close agreement with the findings of Anderson (2002).
Looking more closely at the results shown in Fig. 4, the relief effectively decays during the first 5-6 climate cycles (0.5-0.6 Myr).The average weathering rate also generally winds down and sediment transport gradually increases, although both are affected by temperature variations (Fig. 4e).
After 0.6 Myr, weathering rates and sediment transport vary in a cyclic manner that mimics the change in temperature.The rates of sediment transport and frost cracking are, however, out-of-phase (Fig. 4e), as the cracking rate peaks during cold periods when the MAT maximises FCI, and sediment transport is most efficient during the warmest periods.This is a result of the temperature offset between the most efficient regimes of the two processes, which is ∼ −5 • C for frost cracking and ∼ 0 • C for frost creep (Fig. 3).Creep is thus most active when MAT approaches 0 • C, whereas frost cracking is more efficient at lower temperatures.The achievement of a quasi steady state in the last stages of the simulation is signalled by an almost uniform sediment thickness (2-3 m) across the surface (Fig. 5).This outcome reflects the details of how sediment production depends on sediment thickness, as the exposed bedrock along the grid boundaries and the sediment-covered bedrock in the grid interior now sees more or less the same rate of frost cracking.A stable sediment cover is hence promoted by the humped shape of the relation between sediment production and sediment thickness because this shape allows the same rate of Earth Surf.Dynam., 3, 463-482, 2015 www.earth-surf-dynam.net/3/463/2015/frost cracking to occur for two different sediment thicknesses (Anderson, 2002;Andersen et al., 2015).We note that steady state is not fully attained because the temperature is forced to vary over time.

The influence of temperature
To investigate the influence of temperature, we repeated the simulation described in Sect.3.1, but with climatic cycles varying within a colder temperature interval, −12 ≤ MAT ≤ −6 • C, and a warmer temperature interval, 0 ≤ MAT ≤ 6 • C, respectively.As with the first simulation, a smooth convex surface develops after a few million years in both cases (Fig. 6), but compared to the initial (intermediate) model (−6 ≤ MAT ≤ 0 • C, Fig. 6b), the colder (Fig. 6a) and the warmer (Fig. 6c) scenarios erode more slowly and produce a higher, flatter surface with less curvature.In the cold case, surface erosion is clearly limited by inefficient sediment transport: average temperature is too cold for optimal sediment transport or frost cracking, and so rates of both processes peak during the warmest interglacial periods (Fig. 6a).
In the warmest scenario, frost cracking occurs mainly during the coldest intervals.
The temperature dependence of frost cracking and frost creep causes the sediment cover thickness in each of these three models to undergo repeated cyclic variations in response to the climate forcing (Fig. 7).For the coldest scenario (−12 ≤ MAT ≤ −6 • C), the temperature, frost cracking, and frost creep are all in-phase because both processes are most active when the climate is warmest.Owing to the controls on frost cracking, the sediment cover thickens slightly when temperature rises and thins when climate cools (Fig. 7).However, the mean sediment thickness remains relatively high because sediment creep is generally inefficient, whereas frost cracking is still active, albeit slow, under several metres of sediment cover (Fig. 3a).
The sediment cover is much thinner for the positive MAT situation (0 ≤ MAT ≤ 6 • C) and variations through time are < 0.2 m (Fig. 7).In this case, sediment production is slow because frost cracking is impeded once a thin sediment cover has formed (Fig. 3).Still, the frost-cracking rate increases slightly when MAT approaches 1-2 • C, causing the sediment cover to grow a little during the coldest periods.Sediment creep is almost steady and independent of MAT owing to the comparatively thin cover (Figs.3b and 6c).
In contrast to the colder and warmer models, the average sediment thickness varies significantly for the intermediate model (−6 ≤ MAT ≤ 0 • C).This is a result of the MAT fluctuation between those temperatures optimal for frost cracking and those optimal for frost creep.The sediment cover thickens when MAT drops below −3 • C and thins when frost creep dominates for MAT > −3 • C (Fig. 7).

Weathering-limited versus transport-limited erosion
In the simulations presented so far, rates of surface erosion vary between 0 and 100 m Myr −1 , with long-term average rates of around 10-20 m Myr −1 (Fig. 6).Two processes dictate these rates: (1) the rate of bedrock weathering by frost cracking and ( 2) the rate at which sediment is transported down-slope.The latter is important because the rate of frost cracking depends on sediment thickness (Fig. 3a).
Owing to a fuller understanding of the physics involved, rates of frost creep are, however, better constrained than rates of frost cracking.As documented by Andersen et al. (2015), creep rates are calculated from the modelled variations in water fraction causing the expansion and contraction of sediment during freeze-thaw events due to the growth of segregation ice lenses, and rates are scaled by measurable parameters, such as thermal properties and the sediment expansion ratio, β (Andersen et al., 2015).On the other hand, rates of frost-driven weathering are scaled with the frost-cracking intensity integrated throughout a year, but the scaling parameter k e is not known.It seems intuitive that k e should at least vary between bedrock lithologies, which may affect weathering rates through differences in, for example, tensile strength and fracture density.
To investigate the influence of k e on the modelled weathering rates, we repeat the first experiment (−6 < MAT < 0 • C) several times with different values of k e , ranging over 2 orders of magnitude from 10 −4 to 10 −2 • C −1 a −1 .For each experiment we record the average weathering rate over the 4 Myr simulation period.For low values of k e , the modelled erosion is clearly limited by the rate of frost cracking (weathering-or production-limited), and average erosion rates increase with k e (Fig. 8).However, for higher values of k e (> 0.003 • C −1 a −1 ) the average erosion rate levels out and an increase in k e does not lead to faster erosion.At this point the erosion rate is limited by the rate of frost creep (transportlimited), which in turn is controlled by surface slope and κ (Fig. 3, Eq. 4).The curvature and average slope of the surface increase with k e (Fig. 8) because higher levels of sediment flux are needed to evacuate the sediments produced by relatively rapid frost cracking.The topographic slope therefore increases in order to deliver faster evacuation of the sediment.
To further explore the transport-limited condition (the plateau in erosion rates obtained for high k e values), we repeated the k e sensitivity analysis using two higher values (4 and 6 km) for the grid width L (initially 2 km).Increasing the grid width widens the landform modelled and decreases the average surface slope.As expected, the reduced surface slope decreases the sediment flux and lowers the transport-limited erosion rates (Fig. 8).
Bedrock outcrops are persistent features of the weathering-limited condition (low k e values).The outcrops remain along the central ridge, as tors, or along the edges of the surface where the sediment cover thins to 0. The outcrops are stable, steady-state forms that, once established, erode equally as fast as the surrounding landscape.The main difference concerns sediment thickness.While bedrock outcrops are free of sediment, the surrounds are covered to a depth of several metres, creating a highly bimodal sediment distribution (Fig. 9).The frost-cracking rate is uniform in this setting due to the relation between frost cracking and sediment thickness taking the form of a humped function for MATs ≤ 0 • C (Fig. 3a).As noted above, the humped function predicts frost cracking to be most effective under a finite thickness of sediment, with lower cracking rates accompanying both thinner and thicker sediment covers.Two different sediment thicknesses can hence lead to the same rate of frost cracking.

Experiment 2: a high-relief landscape
The second experiment illustrates the influence of frost cracking and frost creep in a larger landscape with kilometrescale relief.Compared to the previous experiment, we now use a more realistic climatic forcing resembling late Cenozoic cooling (Fig. 10a).The climate forcing is constructed via linear transformation of a marine δ 18 O record (Zachos et al., 2001) to sea-level temperature.So as to calibrate the temperature curve over the last 14 Myr and to approximate the first-order climate history of the North Atlantic margins, we set the present-day sea-level MAT off southern Norway In contrast, the sediment distribution is highly bimodal in the production-limited case.There is no sediment along the edges and on tors that protrude up through the sediment cover along the central ridge.In the other parts of the landscape sediment covers are mostly 2-3 m thick.
coupled with that estimated for the last glacial maximum.The resulting temperature history is inevitably a very rough estimation of past climate; still, it enables us to incorporate the gradual cooling as well as the increasing variability in the climate throughout the late Cenozoic.
Again, we use a simulated fluvial-style initial topography with a central ridge up to 2 km high and fluvial valley structures; this time on a 25 km × 50 km grid where the left and right boundaries are at sea level (i.e. 100 × 200 cells, with a spatial resolution of 250 m).The inclusion of regional isostasy incorporates the long-term effects of unloading the lithosphere by erosion.The isostatic displacements are assumed to be uniform across the grid, which, due to the flexural rigidity of the lithosphere, is a reasonable assumption for the spatial scales used here.Isostasy therefore simply raises the topography uniformly in response to the average erosion of the landscape: where ρ r = 2900 kg m −3 is bedrock density, ρ s = 2400 kg m −3 is sediment density, and ρ a = 3250 kg m −3 is the density of the asthenosphere materials that provide the isostatic compensation.ē is the bedrock erosion averaged across the grid, and S is the change in average sediment thickness.The sediment is able to escape the landscape across the left and right grid boundaries.We set k e = 10 −3 • C −1 a −1 (i.e. the mid-range of k e values applied in experiment 1).The MAT is assumed to decrease linearly with elevation at a lapse rate of 6 • C km −1 , causing the highest peaks to be about 12 • C colder than the lowest parts of the landscape.Because of these internal temperature contrasts, the various parts of the landscape differ in response to the gradual cooling over the model simulation.The frost processes first attack the highest summits where temperatures became sufficiently low for frost cracking to take hold more than 10 Myr ago.Unlike the previous simulations, the initial relief is much too high to be completely erased by the frost-driven processes.However, the combined action of frost cracking and frost creep rapidly transforms the originally sharp summit crests into smooth, high-elevation summit flats in the upper parts of the landscape (Fig. 10b-e).In this experiment, the smoothing of the highest summits takes place primarily before the pronounced cooling in the Quaternary, but the same smoothing effects then spread to lower elevation summits after about 3 Myr ago.
The average erosion rate across the landscape increases throughout the late Cenozoic simulation period, from an initial ∼ 2 to about 6 m Myr −1 , which represents a Pleistocene average, and reaching ∼ 10 m Myr −1 during the coldest periods (Fig. 11a).Within this overall trend, however, we note marked differences between the high and the low parts of the landscape.At low elevations the erosion rate increases dramatically in the latest and coldest part of the Cenozoic (red curve in Fig. 11c), whereas the high-elevation summits experience stagnating or even decreasing erosion rates in the late Quaternary (green curve in Fig. 11c).
Frost-driven processes erode up to 300 m from the highest peaks over the 14 Myr simulation period, whereas ∼ 100-200 m is removed from the lower summits around 1000 m above sea level (Fig. 10).Slope inclinations on and near summits generally decrease in line with the process of flattening at high elevations.The total isostatic rock uplift amounts to 32 m; thus, only areas that erode less than 32 m experience a net surface uplift due to isostasy.

Experiment 3: the influence of alpine glaciations
In the final experiment, we repeat the previous simulation, but introduce alpine glaciations over the last 3 Myr of the simulation.Our goal is to explore how major modifications of the landscape by glacial erosion influence the overall efficiency of frost cracking and frost creep.To form the ice, we apply a simple mass-balance function that relates the rate of ice accumulation and/or ablation, ṁ, to temperature: where T ELA = −2 • C is the temperature at the equilibrium line altitude (where ablation balances accumulation); α = 0.5 m a −1 • C −1 is the accumulation gradient; and β = 1.5 m a −1 • C −1 is the ablation gradient.As before, MAT decreases linearly with elevation at a lapse rate of 6 • C km −1 .After ṁ is calculated for each cell, the ice accumulation, ṁ > 0, is modified to account for avalanches and snow drift by enforcing ice-accumulation penalty functions, ∇ 2 h and ∇h , for topographic curvature and slope, respectively: where ensures that total ice accumulation is maintained during this operation.That is, ice removed from steep or convex areas accumulates elsewhere.
We use the iSOSIA (integrated second-order shallow ice approximation) ice model (Egholm et al., 2011(Egholm et al., , 2012) ) to compute the flow of ice over the final 3 Myr of the simulation.iSOSIA is a depth-integrated model that computes depth-averaged horizontal velocities and rates of basal sliding.The depth-integration aids computational efficiency and enables us to model ice flow over 3 million years with time steps down to a few days.iSOSIA includes longitudinal and transverse stress components, which grant improved accuracy over standard shallow-ice approximations in steep and rugged terrain (where ice-flow velocity varies over short distances) but require iterative loops to solve the non-linear relations between stress and velocity.Subglacial erosion rate is assumed to scale with the rate of basal sliding, which occurs only where the temperature at the ice bed is at the pressure melting point.Consequently, subglacial erosion does not occur where the ice is cold-based.The subglacial thermal state is largely a function of surface air temperature and therefore shifts through time reflecting variations in climate.The reader is referred to Egholm et al. (2011Egholm et al. ( , 2012) ) for further details on the ice-flow model.What is important to note here is that the ice model is designed to produce the following behaviour, which we consider as model input: 1.The extent of the glaciers varies through time according to the climatic forcing.
2. Subglacial erosion lowers the valley floors where the ice is warm-based, and this increases the relief, steepens the valley sides, and accelerates isostatic uplift.

Subglacial erosion and sediment transport do not occur
where the ice is cold-based.
4. Subglacial sediment is transported wherever ice is sliding, which effectively removes sediment covers in glaciated valleys over timescales of 10 4 years.
5. An ice cover thicker than 20 m halts frost cracking altogether because daily and annual temperature variations in the bedrock are dampened by the overlying ice.
The simulated landscape is significantly modified by the inclusion of glaciers (Fig. 12).Almost the entire landscape is covered by ice masses during the coldest glaciations, whereas warmer interglacials are without ice or support only small cirque-type glaciers high in the landscape.Deep and relatively wide glacial troughs develop in the lower parts of the largest catchments; the troughs are over-deepened down to 500 m below sea level.The smaller catchments become hanging valleys.
The summit flats are notably unaffected by subglacial erosion.Where they occur close to valley troughs, some lowering and steepening is evident, but the summit flats and their sediment covers are otherwise stable.This behaviour is not surprising as it stems directly from the prescribed accumulation patterns that minimise ice growth on convex summit flats.When ice does grow over summit flats, it is cold-based and therefore non-erosive.
It is thus more interesting to note the indirect consequences of the subglacial erosion for the rates of frost cracking.Glaciers influence the conditions for frost cracking in several ways: 1.The presence of glaciers stalls frost cracking because daily and annual temperature variations are strongly dampened by the ice cover.
2. Warm-based glaciers work to strip the landscape of sediment, and this may increase frost cracking upon deglaciation where frost cracking would otherwise be transport-limited.
3. Subglacial erosion steepens the valley sides, and the higher surface slopes boost sediment transport rates and therefore also transport-limited frost cracking.
4. Although focussed in the valleys, the average subglacial erosion is substantial and the isostatic response amounts to more than 150 m in this experiment.The slow shift in elevation causes less frost cracking at the highest summit flats where MATs are already lower than those optimal for efficient frost cracking, whereas frost cracking is enhanced at lower elevations due to isostatic uplift lowering temperatures.
The positive implications of subglacial erosion for the periglacial processes are by far the strongest in our experiment: average periglacial erosion rates accelerate markedly in response to the glaciations (Fig. 13a).The glacially steepened valley sides are the first to experience enhanced frost cracking, shortly after deglaciation when freshly scoured glacial valleys are without sediment cover (Fig. 13b).Sediment transport is effective along the steepened valley sides, and the valleys' rates of frost cracking are persistently high during the interglacial periods.In comparison with experiment 2, we note that the inclusion of mountain glaciers results in intensified periglacial erosion along steepened valley sides (up to 400 m) and slightly less (∼ 30 m) bedrock erosion of the highest summit flats (Fig. 13d), which now experience a ∼100 m net rise in surface elevation due to isostasy.

Discussion
The computational experiments presented here integrate rates of periglacial frost cracking (weathering) and frostdriven creep over millions of years and predict landscape evolution based on those processes.The weathering rates are estimated from thermal models that quantify frost-cracking intensity in bedrock based on temperatures, thermal gradients, and the availability of liquid water (Walder and Hallet, 1985;Hales and Roering, 2007;Matsuoka, 2008;Anderson et al., 2012;Andersen et al., 2015).Likewise, the thermal models are used to predict rates of frost creep by integrating the frost-heave activity in sediments on annual timescales (Andersen et al., 2015).As demonstrated by our three experiments, periglacial activity primarily leads to high-elevation summit flats mantled with a veneer of sediment.These smoothly convex surfaces are emergent properties of periglacial surface processes, as shown by previous studies (Anderson, 2002;Anderson et al., 2012).In some of our model scenarios the rise of bedrock outcrops (tors) denotes weathering-limited conditions where sediment transport matches, or exceeds, the tempo of bedrock weathering (Fig. 8).In such settings, a bimodal arrangement develops in which bare bedrock is flanked by sediment mantles 2-3 m thick (Fig. 9).The bimodal sediment distribution is enhanced by the humped weathering function in cold settings (MAT < 0 • C), wherein weathering rates are maximised beneath a ∼ 2 m thick sediment cover and decline with either thinning or thickening (Fig. 3a).
A key aspect of summit-flat dynamics is the coupled nature of sediment production and transport.Where bedrock weathering rates increase (via higher k e values in our model), the system quickly becomes transport-limited as bedrock weathering nearly stalls once the cover exceeds about 3 m thickness (Fig. 3a).This feedback reflects the presence of watersaturated sediment that impedes the penetration of diurnal and/or annual temperature variations into the bedrock; an effect enhanced by latent heat in the sediment.Our results suggest that transport via frost creep, which in our model depends on temperature, sediment thickness, and frost susceptibility (parameter β in Andersen et al. (2015)), limits longterm frost weathering to rates of the order 10 m Myr −1 for landscapes with low to moderate slopes.Such denudation rates are low but in good agreement with those estimated via cosmogenic nuclide measurements on cold-region blockfield and/or summit flats (Small et al., 1997(Small et al., , 1999;;Goodfellow, 2012).In our model, the mode of sediment transport was restricted to frost creep so we could explore its effectiveness in isolation.However, noting the transport-limited nature of summit-flat regolith highlights the inherent conservatism  12, which combines subglacial and periglacial erosion (over the final 3 Myr of the total 14 Myr simulation).The red curve is subglacial erosion rate, and the green curve is the periglacial erosion rate.The black curve shows the periglacial erosion rate from the experiment without glaciers (same as black curve in Fig. 11a).The blue curve is sea-level MAT.The difference between the black and green curve highlights the effect of subglacial erosion on the periglacial erosion rate, which is of our model predictions.The inclusion of other sediment transport mechanisms, such as stream flow and bioturbation, would serve to accelerate summit lowering.The transportlimited nature of the erosion may also explain why frost cracking shows much greater efficiency (rates ∼ 1 mm a −1 ) in steeper, tectonically active landscapes where thick mantles neither develop nor persist in the long term (Hales and Roering, 2009;Delunel et al., 2010).

Assumptions and limitations
The simulated landscape evolution is a direct consequence of the assumptions in our model.We have employed the latest understanding of the physical principles that underpin large-scale periglacial landscape evolution (Hales and Roering, 2007;Anderson et al., 2012); however, many different assumptions could be made, which would perhaps lead to different results.In the following we attempt to outline some of the main limitations of our approach.

Substrate factors
The parametrisation of frost creep depends on the ability of the bulk sediment to expand when freezing and contract when thawing (parameter β in Andersen et al. (2015)).This process is boosted by the presence of densely packed finegrained (silt-sized) sediment with fully water-saturated pore spaces (Chamberlain, 1981).However, it remains open question how effectively and at what rate cold-region weathering processes produce the silt-sized grains.Importantly, the rates of transport and erosion may be reduced if weathering processes do not produce small grains quickly enough.We note that the abundance of fine-grained material observed on many summit flats today (e.g.Strømsøe and Paasche, 2011;Goodfellow, 2012) might indicate that multiple modes of frost-driven sediment production and/or chemical weathering are keeping pace with transport processes under current conditions.However, it is not straightforward to deduce a general case for the long timescales involved in our model simulations.
Rock types clearly differ in their susceptibility to frost cracking (Lautridou and Ozouf, 1982;Hallet et al., 1991;Matsuoka, 2001), and although the model disregards this issue, our results imply that lithology will be important for periglacial landscape evolution (Goodfellow et al., 2014).The potential for water to migrate towards incipient ice lenses first of all requires some degree of bedrock permeability; frost cracking thus depends on pre-existing fractures in addition to bedrock porosity.In summarising the results of frost-cracking experiments with various rock types, Matsuoka (2001) reports frost cracking for high-porosity rocks (tuff, shale, chalk) from −5 to 0 • C and for medium-porosity rocks (limestone, sandstones) between −6 and −3 • C. Lowporosity bedrock, on the other hand, shows little or no cracking even under optimal moisture conditions (Matsuoka, 2001).Together, these observations indicate that frost susceptibility may be highly lithology-dependent.The simplest way to incorporate lithological differences in our model would be to vary k e between grid cells in the simulated landscape.

Water availability
Moisture is often present in mountain environments, but in certain cases precipitation may be more important than temperature in limiting frost action (Sass, 2005;Hall and Thorn, 2011).Hyper-arid regions, for instance, are too dry to promote frost cracking (Hall et al., 2002) and polar deserts sustain the slowest denudation rates on Earth (Portenga and Bierman, 2011).In less extreme environments it is worth noting that water need not be present all year to promote frost cracking (Girard et al., 2013).For negative MAT environments, water must be available at the surface in warm periods (e.g. during spring when surface warming drives snow melt or during summer when precipitation falls as rain).For positive MAT (periglacial) environments, water must be present during winter at depth within the bedrock when the surface temperature drops into the frost-cracking window.This is expected to apply in areas where autumn and/or winter are wet, or where the bedrock is prevented from drying out due to low insolation or local topographic factors.
Quantitative measurements of rock moisture in natural environments are scarce and mainly refer to steep rock faces with little or no sediment cover.Existing data point to quite consistent moisture levels at depths of a few decimetres (Sass, 2005).The observations generally highlight the importance of direct precipitation and local climatic properties, such as prevailing wind direction, insolation patterns, and the distribution of snow, but it is doubtful that such findings apply to regolith-mantled bedrock.Our model experiments assume total water saturation in both bedrock and sediment at all times, without fluctuation.This assumption probably yields an overestimation of frost action.However, we note that an overall reduction in process rates will not significantly alter model outcomes but only affect the timescale of landscape change.On the other hand, prevailing wind direction and insolation patterns can, by influencing temperature and water availability, have aspect-dependent effects that lead to asymmetrical denudation patterns (Anderson et al., 2012).Our model simulations do not include such factors.

The development of summit flats
There are many aspects of periglacial processes that our landscape evolution simulations do not address directly.However, the emergence of flattish or gently convex areas (summit flats) represents a robust and highly reproducible outcome of all our experiments.For example, the emergence of summit flats is independent of inherent assumptions regarding how water availability affects FCI (Fig. 14).In fact, low-relief surfaces are characteristic of all types of transport-limited weathering, irrespective of climate.They follow from diffusive, slope-dependent transport processes that are a function of sediment thickness and the absence of a channel (advective) network.Once sediment cover is present via whichever formation process (frost cracking, chemical weathering, or aeolian accession), the transport dynamics differ little and operate according to the diffusionlike mechanics in the model.Regolith-mantled summit flats therefore seem to arise as an inevitable consequence of longterm transport-limited behaviour on hillslopes.
A distinguishing attribute of the periglacial environment lies in the strong temperature dependence of frost cracking, which leads to weathering and erosion to be especially prevalent at particular elevation intervals.When such processes are coupled with slow rates of denudation, frost action may be concentrated on the same portion of the landscape for a long period.The slow erosion rates on summit flats predicted by our model and corroborated by cosmogenic nuclide measurements suggest that pronounced summit flats develop primarily in settings where other erosional agents, especially rivers and glaciers, are subdued relative to frost action and where cold climatic conditions have prevailed over millionyear timescales.These criteria are likely met in tectonically quiescent regions at mid to high latitudes.
To contextualise the low denudation rates predicted by our model, we note that a recent compilation of apatite fission track thermochronometry suggests that the average rates of post-Devonian erosion in the Scandinavian mountains have been less than 10 m Myr −1 (Medvedev and Hartz, 2015).Such low rates are also in agreement with average millennialscale bedrock erosion rates derived from cosmogenic nuclides in many areas of the planet (Portenga and Bierman, 2011).The range of measured erosion rates confirm that average erosion rates on Earth are indeed low, with clear exceptions in actively deforming orogens (Montgomery and Brandon, 2002;Von Blanckenburg, 2005) and under fast-moving glaciers (Hallet et al., 1996).The correspondence between empirical erosion rate data and the outcomes of our computational experiments, in turn, supports the argument that the slow frost-driven processes simulated here could in many settings represent the primary agents of landscape denudation.

Implications for glaciated passive continental margins
We now consider the implications of the model experiments for the evolution of glaciated passive margin landscapes (e.g. in Norway and Greenland), which bear a resemblance to the modelled scenarios (Figs. 10 and 12).Our results indicate that the highest (coldest) summit flats were largely established prior to deep cooling in the Quaternary (Fig. 10), and such areas have now become too cold for frost-driven processes to operate effectively -even during interglacial periods.The apparently minimal present-day frost action observed on the highest summit flats in Norway (Strømsøe and Paasche, 2011) is in agreement with our model results, which predict the summit flats slowly evolving into landscapes with extremely low rates of geomorphic activity due to global cooling and isostatic surface uplift.Low-relief areas at high elevations have been traditionally interpreted as evidence for accelerated rock uplift following peneplanation at base level (Lidmar-Bergström et al., 2000;Japsen et al., 2009).Our experiments provide a viable alternative: low-relief areas may develop more or less in situ via mechanisms unrelated to either base level or tectonism (Anderson, 2002;Nielsen et al., 2009;Steer et al., 2012).Whether high summit flats in the mountains of Norway and Greenland are uplifted remnants of peneplains (Lidmar-Bergström et al., 2000;Japsen et al., 2009, e.g.) or the product of long-standing topography and global cooling (Nielsen et al., 2009;Steer et al., 2012) needs further investigation.Nonetheless, diffusional, transport-limited processes can undoubtedly produce smoothly convex hillslopes and summits (Anderson, 2002); therefore, low-relief morphology alone is insufficient evidence for invoking tectonism or for identifying remnants of fluvial landscapes formed under "preglacial" conditions (Kleman and Stroeven, 1997;Hall et al., 2013;Lidmar-Bergström et al., 2013).
The highest summit flats in our simulations experienced decelerating frost cracking during the Quaternary mainly due to extremely cold conditions stalling rates of sediment creep (Fig. 11c, green curve).However, at the same time, global cooling over the last ∼ 3 Myr introduced frost action to lower elevations, which represent a larger proportion of the mountain topography.Hence, the landscape as a whole experienced significant intensification of frost activity (Fig. 11a).This trend naturally depends on the specific temperature variations, the relationship between temperature and elevation, as well as the landscape hypsometry.Nevertheless, a major areal expansion of Earth's periglacial realm suggested by our results is notably compatible with the proposed acceleration of sediment production in the late Cenozoic (Hay et al., 1988;Zhang et al., 2001;Molnar, 2004;Herman et al., 2013).Global cooling expanded both glacial and periglacial activity.Figure 14.Summit-flat evolution using FCI patterns resulting from different assumptions regarding water availability (see also Andersen et al., 2015, Fig. 11).Panel (a): frost cracking occurs whenever water is present in the direction of increasing temperature (resembling the model by Hales and Roering, 2007).Panel (b): the distance from site of frost cracking to the nearest water source (T > 0 • C) scales FCI in accordance with the model of Anderson et al. (2012).Panel (c): like (b), but now porosity and water fraction at the water source affect FCI as in Andersen et al. (2015).Compared to our standard model (Fig. 3), (c) has constant flow restriction parameters (2 m −1 ).Apart from the difference in FCI, all three model simulations are identical to our experiment 1 (Sect.3.1).The differences in how the three FCIs depend on sediment thickness clearly influence the sediment thickness.However, although surface curvature and final sediment distribution vary, all three models inevitably result in the development of a summit flat.
we suggest that periodic shift between frost-driven erosion and subglacial erosion was potentially a key driver of the accelerated denudation observed in many mid-and highlatitude mountains.

Conclusions
The landscape evolution models presented here integrate the action of frost cracking and frost creep over million-year timescales.The rates of cracking and creep are based on thermal models that record temperature gradients and frost-heave events in the subsurface (Andersen et al., 2015).We have devised three types of computational experiments that focus particularly on the development of high-elevation summit flats in a periglacial environment on different spatial scales.The third experiment involves both periglacial and glacial erosion processes.The key findings are as follows: 1. Frost cracking and frost creep inevitably lead to smoothing of the relief on summits and hillslopes.This is a highly reproducible and robust outcome based on a wide range of temperatures, provided that atmospheric conditions drive occasional freezing below −3 • C together with intervals of positive temperatures that allow water to be mobilised.
2. Owing to the climate dependence of periglacial activity and atmospheric cooling at higher altitudes, the efficiency of frost action is optimised in certain elevation intervals.This leads to the formation of smooth surfaces at altitudes that have experienced climatic conditions optimal for periglacial processes over long time periods.This result indicates that low-relief areas can form at high elevation due to surface processes and that flatness alone is not diagnostic of a formation history close to sea level.It is important to note that the modelled topographic smoothing takes place on timescales of millions of years, in agreement with erosion histories estimated from such settings in nature (Small et al., www.earth-surf-dynam.net/3/463/2015/Earth Surf.Dynam., 3, 463-482, 2015Dynam., 3, 463-482, 1997;;Fabel et al., 2002;Phillips et al., 2006;Goodfellow et al., 2014).
3. Introducing a simple scaling of temperature to marine δ 18 O records so as to mimic the late Cenozoic climatic cooling cultivates the notion that surfaces shaped by periglacial processes in the Miocene and Pliocene can have experienced declining activity through the Quaternary due to deep cooling.This may potentially explain the low activity observed on some summit flats today (e.g.Rea et al., 1996;Whalley et al., 2004;Strømsøe and Paasche, 2011;Goodfellow, 2012).
4. Models incorporating both periglacial and glacial erosion potentially explain the highly bimodal topography observed in many alpine settings, with glaciers rapidly carving out deep valleys separated by low-relief interfluves strewn with slowly evolving periglacial regolith mantles.
5. Although the high summit flats may have experienced slower frost-driven erosion during the coldest Quaternary intervals, average rates of frost cracking across large areas are likely to have accelerated over the last few million years.The experiment combining glacial and periglacial erosion further suggests that subglacial erosion may have increased periglacial erosion significantly during interglacials by steepening valley sides and removing sediment cover via warm-based subglacial erosion.

Figure 3 .
Figure 3. Frost-cracking intensity (FCI) (a) and frost-driven sediment transport efficiency (κ) (b) contours shown as functions of mean annual temperature (MAT) and sediment cover thickness (S).The contour plots are created from modelled values of frost cracking and frost creep for a range of combinations of MAT and S.These results were obtained by solving the one-dimensional heat equation and integrating the rates through annual temperature cycles.The model is documented in detail in the accompanying paper(Andersen et al., 2015).

Figure 4 .
Figure 4. Computational experiment (no. 1) with periglacial erosion of a synthetic, fluvial-style initial landscape.The experiment shows the development of a single summit flat.The model landscape is exposed to 20 recurrent climate cycles of 100 kyr each, entailing saw-tooth variations in temperature.Panel (a): initial fluvialstyle topography.Panels (b-d): predicted topography after 0.1, 0.5 and 2 Myr, respectively.Bottom panel time series show temperature (blue), average erosion rate (grey), and average frost-creep diffusivity (black).

Figure 5 .
Figure 5. Same model experiment (no. 1) as in Fig. 4 but extended to 4 Myr and showing average relief, topography, and sediment thickness over time.The model evolution is shown at 0, 0.5, 2, and 4 Myr from left to right.Top panel: the topographic envelope of the landscape.Green area marks the difference between the minimum and maximum topography when measured along the long axis (3 km) of the grid (black line indicates mean elevation).Note that the green area collapses over time because the initial local relief is removed by the periglacial processes.Middle panel: modelled topography in map view.Bottom panel: distribution of sediment in the simulated landscape.Note that the initial relief decays relatively rapidly, but traces of the initial relief are still preserved in the sediment distribution after 2 Myr.However, the sediment thickness becomes uniform (3-4 m) in the final stage of the simulation (4 Myr), and the modelled landscape no longer carries any memory of its initial configuration.

Figure 6 .
Figure 6.Development of a summit flat under three different climatic scenarios: (a) a cold model (−12 < MAT < −6 • C); (b) an intermediate model (−6 < MAT < 0 • C, same as Fig. 5); and (c) a warm model (0 < MAT < 6 • C).Final topography at 2 Myr is shown in left panels.Right panels show time series of temperature (blue), average frost-cracking rate ( ė, grey), and sediment diffusivity (κ, black) over four climate cycles between 1.6 and 2.0 Myr after simulation start.Note that frost cracking and frost creep are in-phase for the cold model but out-of-phase for the intermediate and warm models.

Figure 7 .
Figure 7. Evolution of mean sediment thickness (red curves) on the three summit flats shown in Fig. 6.Rates of frost cracking as function of MAT and sediment thickness are shown as background contours.Letters a, b, and c refer to the summit flats shown in Fig. 6.Note that in all three cases mean sediment thickness varies in a cyclic manner when MAT increases and decreases over each climate cycle.The cold and the intermediate models have much thicker sediments (∼ 3 m) than the warm model (∼ 0.5 m).The intermediate model shows the strongest thickening and thinning of the sediment cover because the MAT varies between the two zones of most effective frost cracking and most effective frost creep.

Figure 8 .
Figure 8. Mean accumulated frost-driven erosion after 4 Myr for different values of k e (frost-cracking scaling factor).Results shown are for the intermediate climate model with −6 < MAT < 0 • C. The three curves are for model grids of different widths.The landscape is steepest for small L. Frost cracking increases with k e when erosion is production-limited.The three curves reach asymptotic levels for larger k e values because erosion becomes transport-limited.The insets show some of the modelled summit flats after 4 Myr.The curvature and slope of the summit flats depend on the rates of erosion because slopes contribute to the scaling of sediment transport.Tors and bedrock ridges are stable phenomena in the left-most landscape where erosion is slow and production-limited.

Figure 9 .
Figure 9. Results for a transport-limited (top) and a production-limited model (bottom).k e = 0.005 • C a −1 for the transport-limited model and 0.0001 • C a −1 for the production-limited model.This is the only difference between the two model experiments, which both have −6 < MAT < 0 • C and L = 2 km.Left panels: modelled topography after 4 Myr.The contour spacing is 10 m.Middle panels: sediment thickness in the landscape after 4 Myr.Right panels: histogram of sediment thickness for the two cases.Almost all of the landscape has > 2 m sediment in the transport-limited case, and sediment thickness gently decreases from 3.5 m in the centre to about 2 m along the edges.In contrast, the sediment distribution is highly bimodal in the production-limited case.There is no sediment along the edges and on tors that protrude up through the sediment cover along the central ridge.In the other parts of the landscape sediment covers are mostly 2-3 m thick.

Figure 10 .
Figure 10.Computational experiment (no.2) based on a landscape with kilometre-scale relief.Panel (a): the 14 Myr temperature curve used as input for the model (left), and the initial topography (right).Panels (b-e): topographic evolution (left panel) and elevation-slope relations (right panel) after 0, 4, 10, and 14 Myr, respectively.To magnify the detail, topography is only shown for a selected region, as outlined by the black rectangle in the right panel of (a).The elevation-slope diagrams each have two sections: upper section shows the distribution of surface slope as box-and-whisker plots for different elevation intervals; lower section shows the frequency of flat regions (slope < 0.1) at the different elevations.Over time, the highest parts of the landscape are smoothed significantly and summit flats are formed at several levels above 1000 m.The formation of summit flats is accompanied by a distinct decrease in surface slope in the highest part of the landscape, and above 1500 m almost 100 % of the landscape has surface slopes less than 0.1.The smoothing starts at elevations above 1500 m (c), but spreads to elevations above 1000 m as the climate cools further during the last interval mimicking the Quaternary (d, e).

Figure 11 .
Figure 11.Long-term erosion outcomes for experiment (no.2) shown in Fig. 10.Panel (a): sea-level temperature (T sl , blue) and average erosion rate ( ė, grey) spanning the 14 Myr simulation period.Erosion rates are shown to vary intensely with temperature; the overall trend is highlighted by a running mean (solid black line).Average erosion rate decays through the first 14 Myr in response to the gradual smoothing of the summit flats.The effects of cooling takes over in the final 6 Myr when average erosion rates accelerate markedly as frost action spreads to lower elevations.Panel (b): total erosion after 14 Myr.Panel (c): erosion rates measured in the two points marked by arrows in (d).Both points are at local summits but differ in elevation and therefore experience different temperatures.The high summit (green curve) has relatively high erosion rates at first, but erosion slowly stalls in the final 6 Myr because the climate becomes too cold for efficient sediment transport by frost creep.In contrast, the lower summit (red) experiences a significant acceleration of erosion in the Quaternary.

Figure 12 .
Figure 12.Computational experiment (no. 3) combining periglacial and subglacial erosion.The boundary conditions are the same as in Fig. 10, with the addition of glaciations over the last 3 Myr of the total 14 Myr simulation.Panel (a): initial fluvial-style topography.Panel (b): periglacial landscape after 11 Myr, just before glaciations start.Panel (c): a snapshot of a period with glaciers (after 13.2 Myr).Note that the ice cover varies greatly with temperature and the final 3 Myr are dominated by glacial cycles that repeatedly cover almost the entire landscape.During the interglacials the landscape is either free of ice or includes small valley glaciers high up.Panel (d): final landscape after 14 Myr.Glaciers have carved out deep troughs that are over-deepened to almost 500 m below sea level.Smaller valleys have become hanging valleys.Summit flats are preserved owing to the presence of non-erosive, cold-based ice on the interfluves.The final topography is highly bimodal, with (1) a high plateau region comprised of many summit flats at similar elevation and (2) deep glacial valleys incised to or below sea level.

Figure 13 .
Figure 13.Panel (a): long-term erosion outcomes for experiment (no. 3) shown in Fig.12, which combines subglacial and periglacial erosion (over the final 3 Myr of the total 14 Myr simulation).The red curve is subglacial erosion rate, and the green curve is the periglacial erosion rate.The black curve shows the periglacial erosion rate from the experiment without glaciers (same as black curve in Fig.11a).The blue curve is sea-level MAT.The difference between the black and green curve highlights the effect of subglacial erosion on the periglacial erosion rate, which is significantly accelerated by the glacial landscape modifications.Panel (b): details of erosion rates during the final 1 Myr, showing subglacial erosion rate (red line), periglacial erosion rate (green line), and sea-level MAT (blue line).Panel (c): final stage of the landscape (same as Fig.12d) with detail magnified at location shown by the black rectangle in Fig.10a.Panel (d): total periglacial erosion.Note that the glacially steepened valley sides are most effectively attacked by frost cracking.
Figure 13.Panel (a): long-term erosion outcomes for experiment (no. 3) shown in Fig.12, which combines subglacial and periglacial erosion (over the final 3 Myr of the total 14 Myr simulation).The red curve is subglacial erosion rate, and the green curve is the periglacial erosion rate.The black curve shows the periglacial erosion rate from the experiment without glaciers (same as black curve in Fig.11a).The blue curve is sea-level MAT.The difference between the black and green curve highlights the effect of subglacial erosion on the periglacial erosion rate, which is significantly accelerated by the glacial landscape modifications.Panel (b): details of erosion rates during the final 1 Myr, showing subglacial erosion rate (red line), periglacial erosion rate (green line), and sea-level MAT (blue line).Panel (c): final stage of the landscape (same as Fig.12d) with detail magnified at location shown by the black rectangle in Fig.10a.Panel (d): total periglacial erosion.Note that the glacially steepened valley sides are most effectively attacked by frost cracking.