Ice flow models and glacial erosion over multiple glacial–interglacial cycles

Mountain topography is constructed through a variety of interacting processes. Over glaciological timescales, even simple representations of glacial-flow physics can reproduce many of the distinctive features formed through glacial erosion. However, detailed comparisons at orogen time and length scales hold potential for quantifying the influence of glacial physics in landscape evolution models. We present a comparison using two different numerical models for glacial flow over single and multiple glaciations, within a modified version of the ICE-Cascade landscape evolution model. This model calculates not only glaciological processes but also hillslope and fluvial erosion and sediment transport, isostasy, and temporally and spatially variable orographic precipitation. We compare the predicted erosion patterns using a modified SIA as well as a nested, 3-D Stokes flow model calculated using COMSOL Multiphysics. Both glacial-flow models predict different patterns in time-averaged erosion rates. However, these results are sensitive to the climate and the ice temperature. For warmer climates with more sliding, the higher-order model yields erosion rates that vary spatially and by almost an order of magnitude from those of the SIA model. As the erosion influences the basal topography and the ice deformation affects the ice thickness and extent, the higher-order glacial model can lead to variations in total ice-covered area that are greater than 30 % those of the SIA model, again with larger differences for temperate ice. Over multiple glaciations and long timescales, these results suggest that higher-order glacial physics should be considered, particularly in temperate, mountainous settings.


Introduction
Over geological time, mountainous topography is formed through a combination of erosional and tectonic processes.In many regions, mountain topography has felt the effects of glacier erosion, in addition to other geomorphic processes.Quantifying the effects of glaciation on topography requires consideration of the physics and rheological properties governing glacial erosion.This study builds upon previous work and evaluates how different assumptions and levels of complexity used in glacial-flow models impact the magnitude of erosion over multiple glacial-interglacial cycles and geologic timescales.This type of study is important for evaluating what level of model complexity (and computational sophis-tication) is required to quantify glacial erosion processes and sediment production in mountainous regions.
Numerical models have often been used to study the influence of glacial erosion on landscape development.These models have ranged from simple 2-D glacial profile models (Anderson et al., 2006;MacGregor et al., 2009) to more complex 3-D models that incorporate a variety of processes and mechanisms (Kessler et al., 2008;Egholm et al., 2012a;Pedersen and Egholm, 2013).Other studies have incorporated glacial erosion into landscape dynamic models that also include fluvial and hillslope erosional processes (Herman and Braun, 2008;Egholm et al., 2011;Yanites and Ehlers, 2012).The evolution and continued development of glacial-flow and erosion models has resulted in the simulation of increasingly Published by Copernicus Publications on behalf of the European Geosciences Union.complex processes such as the influence of subglacial hydrology (Egholm et al., 2011;Herman et al., 2011;Iverson, 2012).Despite these advances, other mechanisms are still represented by simplified assumptions and approximations, particularly the underlying physics of ice flow.Within the field of glaciology, as computing power has increased, higher-order glacial-flow models (Pattyn et al., 2008) have been made more accessible.These higher-order processes are often simulated over the timescales useful for glacial and climatic studies (10 3 to 10 4 years), yet still shorter than the timescale of orogen topographic development and Quaternary glaciations (10 5 to 10 7 years).The incorporation of these higher-order models into orogenic timescale models can be useful to better represent the glacial flow in alpine settings, where the effects of longitudinal and lateral stresses on glacial flow and erosion should be important (Egholm et al., 2012a, b).
In many orogenic-scale models, glacial flow and erosion have been represented using simplifying assumptions, such as the shallow ice approximation (SIA) for glacial flow (Kessler et al., 2008;Iverson, 2012).This approximation simplifies the ice flow equation (Glen flow law) by only considering the first-order simple shear stresses (Cuffey and Paterson, 2010).While this approximation is appropriate for some specific glacial settings, particularly ice sheets where surface slopes are shallow, ice thickness is small compared to ice extent, and sliding velocities are small compared to deformational velocities; for alpine glaciers this assumption misses effects that result from glacial flow through narrow and steep topography (Egholm et al., 2011).While even the simplest approach has its merits, defining the conditions under which a higher-order model should be used warrants more detailed consideration.Recent work has shown that, over the length and timescales of glacial valley formation, higher-order (HO) glacial-flow models have important feedbacks (Egholm et al., 2009(Egholm et al., , 2012b;;Pedersen and Egholm, 2013).Specifically, Egholm et al. (2012a, b) showed that, on the glacier valley scale and over a single glacial cycle, the incorporation of lateral and longitudinal stresses can provide an important mechanism for suppressing potential runaway problems that can come from using the SIA.In addition, larger, regional models have been investigated, and other investigators have stressed that the physics and form of the glacier are certainly important, and that models of alpine glaciers and their erosional patterns are influenced by the choice of physics, particularly as the landforms and valley profiles evolve (Pedersen and Egholm, 2013).
Previous work has addressed differences between SIA and higher-order models, which includes work that focuses on glaciers in steep terrain (Egholm et al., 2011) and glaciers in landscape evolution models (Egholm et al., 2012b).For glaciers in alpine terrain and in landscape evolution scales, the SIA approach was found to be less accurate in predicting sliding velocities and patterns of basal shear stress, since this approximation fails to incorporate lateral and longitudi-nal stresses in glacier valleys.The SIA led to positive feedbacks where enhanced erosion caused deeper and steeper topography that in turn led to higher sliding and erosion rates.hese feedbacks tend to be dampened when used within full Stokes and higher-order models which incorporate at least approximations of the lateral and longitudinal stresses.Other studies have also addressed comparisons between SIA and HO ice flow models, including the benchmarks developed and tested in Pattyn et al. (2008).However, these benchmarks pertain mostly to 3-D continental ice sheets.The benchmarks themselves only include 2-D flow-line models for an alpine glacier and focus only on non-evolving ice flow and basal shear stress.
This study complements previous work by evaluating the effect of the glacial-flow physics model on predicted variations in glacial erosion over both single and multiple glacial cycles.This study can be differentiated from other work due to the following reasons.(1) The timescale investigated in this study is significantly longer, and extends to geologic timescales of 400 000 years, including three full glaciations, with a range of different climate scenarios.Thus, we report here the effects of differences between the two approaches over temporally varying, and multiple, glacial-interglacial cycles.This effect is potentially important because catchment hypsometry evolves over long timescales, as does the thermal state of the glacier.Furthermore, (2) previous work by Egholm et al. (2011Egholm et al. ( , 2012b) used a 3-D second-order shallow ice approach (see also Pedersen and Egholm, 2013) for comparison to a traditional SIA model.A full-stress Stokes flow model was used only in 2-D for comparison to the SIA and could not evaluate 3-D topographic effects on ice mechanics.We evaluate the strengths and weaknesses of SIA and Stokes flow within a 3-D landscape evolution model and address how a HO glacial model might affect topography and its evolution over multiple glacial-interglacial periods.Our end goal is to add to the understanding of when and under what conditions more simplified models, such as the SIA are sufficient, as applicable to larger-scale problems such as sediment production, ice-sheet stability, and tectonic-erosional interactions under alpine glaciers and continental ice sheets.While many simplifications are used in this comparison, such as the choice of erosion and sliding laws, which are discussed more thoroughly in Sect.4.4, this comparison ultimately yields some quantitative evidence of how HO glacialflow physics can influence the evolution of topography.This comparison also provides more evidence that the choice of glacial-flow physics in landscape evolution models should be made at least partially based upon the climate of the model and the timescales of interest, along with the topography and glaciation style (i.e., ice sheets versus alpine glaciers).

Methods
Here we build upon previous work and investigate the influence of glacial-flow physics on a developing orogen over geologic timescales using a modified version of the ICE-Cascade landscape evolution model.In order to compare the importance of the choice of ice physics, simulations are run using both the SIA model and the nested HO model.We repeat these comparisons over different climate scenarios in order to highlight temperature-dependent effects.The different climate simulations are used because glacial erosion is dependent upon the existence of liquid water at the base of the glacier, a temperature-dependent property.For brevity, a summary is given of ICE-Cascade, its major components, the physics governing the ice flow, and the modeling framework behind the HO glacial-flow model.All relevant model parameters used for ICE-Cascade are presented in Table 1; readers are referred to the associated references for additional details.

Simulations
Three separate simulations based on climatic and ice temperature conditions are performed using both the SIA and the HO physics models described in Sect.2.3.These simulations are summarized in Table 2, and the sea-level temperature patterns over time are shown in Fig. 3.
-Simulation 1 uses a sinusoidal temperature pattern with amplitude of 6 • C and sea-level minimum temperature of 2 • C, with a wavelength of 100 kyr.
-Simulation 2 has a similar pattern but with the sea-level minimum shifted to 0 • C, so there are more instances of cold ice where the base is frozen.
-Simulation 3 uses the same temperature pattern as Simulation 1; however, in this simulation sliding occurs everywhere the ice thickness is greater than 10 m, and the ice temperature at the base of the glacier is not factored into this calculation.
Each simulation was run for over 400 kyr, with glaciations major occurring every 100 kyr.Over three full glaciations are captured during this time interval.Figure 1 shows Simulation 2 (SIA) with hillshade topography and the ice coverage and thickness at T = 100 kyr in the simulation.Compared to any individual Pleistocene glacial pulse, the simulated glaciations are extensive and long-lived, particularly with the sinusoidal climate forcing.This is to emphasize long-term evolution effects versus influence from a single, quick glaciation.

ICE-Cascade orogen development model and climate parameters
ICE-Cascade allows us to model topographic evolution over geologic timescales, with the influences of both constructive (tectonics and sediment deposition) and destructive (erosion) processes (Herman and Braun, 2008;Yanites and Ehlers, 2012).At each time step, the topography is uplifted according to an input rate (Table 1) along with a component based on flexural isostasy.The isostatic response of the landscape is affected by loading and unloading due to erosion, sediment deposition, and the ice thickness.Following the uplift, the landscape is eroded according to the rates from the fluvial, hillslope, and glacial modules.Where glacier ice is nonexistent or thinner than 10 m, fluvial and hillslope processes erode and redistribute sediment (Yanites and Ehlers, 2012).Sediment transport by the glaciers occurs in regions of ice thicker than 10 m; bed material is eroded and immediately transferred to the fluvial system that emanates from the toe of the glacier.River discharge is calculated based upon the upstream precipitation and water-equivalent ice melt from upstream regions.Fluvial erosion processes are calculated from this discharge and the sediment supply, local topographic slope, and the channel width, all of which also are input into a linear sediment cover model (Braun and Sambridge, 1997).Hillslope processes are simulated from diffusion and a threshold landsliding algorithm when hillslopes steepen beyond a certain threshold (Burbank et al., 1996;Stolar et al., 2007).Within ICE-Cascade, the climate simulation is a combination of the inputs that govern the pattern of sea-level temperature and an orographic precipitation model (Yanites and Ehlers, 2012;Roe et al., 2003).The temperature and moisture content variations over all elevations is calculated using an input lapse rate and the Clausius-Clapeyron relation, and these values, along with inputs governing the wind speed and direction, are then used to calculate annual precipitation (Roe et al., 2003).When the temperature is below freezing, the precipitation takes the form of snow.A positive-degreeday algorithm is used to determine the number of days above freezing in any given year (Braithwaite, 1995), and this, in turn, is used in calculate the amount of melt of snow and ice.For any point in the landscape, the annual mass balance is then simply the difference between the amount of snowfall and the amount of melt.Climate parameters governing these processes are given in Table 1.
For these simulations, the initial topography is an equilibrium landscape generated using only fluvial and hillslope processes.This topography was built from an earlier simulation (simulation m01 from Yanites and Ehlers, 2012) that started with random topography and was allowed to develop over 16 Myr with the same hillslope, fluvial, uplift, and climate parameters as found in Table 1.

Glacial models
At the beginning of the simulations, glaciations evolve where a persistent positive mass balance exists.Glaciers flow from the orogen and its valleys onto the continental slope, where they form piedmont lobes.The shelf can be seen in Fig. 1, extending to the edges of the y domain from 0 to 70 km and from 225 to 245 km.The shelves, with a slope of 0.001, were added to ensure that the ice velocities at the orogen-parallel boundaries are numerically stable (Yanites and Ehlers, 2012).The shelf edges (y = 0 and 245 km) have Dirichlet boundary conditions, where their elevation is fixed to sea level and does not change throughout the simulations.During these simulations, the ice never reached these boundaries.
The depth-averaged and sliding velocities are calculated from the glacial-flow physics (see Sects.2.3.1 and 2.3.2 for more details).The sliding velocity is only calculated where the temperature at the base of the glacier is at the pressure-melting point associated with the local ice thickness (except in Simulation 3, where there is sliding everywhere).The basal temperature is determined from a conductive heat model for ice, where the upper boundary is the surface temperature, and the lower boundary is the geothermal heat flux (Table 1).Beyond glacial flow, other ice feedbacks are incorporated in the model.On steep slopes avalanching occurs when the snow surface slope exceeds a critical stability angle (Table 1), in which case all the snow at the point is redistributed to the next downstream node.Where the glacier terminates below sea level on the shelf, iceberg calving occurs and is a function of the depth of ice below sea level.In addition to calving, a buoyancy feedback is also incorporated for when the glacier erodes an overdeepening below sea level (Yanites and Ehlers, 2012).In an overdeepened, below-sea-level reach, sea level is assumed to be the equipotential ground water surface, such that the glacial erosion rate decreases as the ice approaches buoyancy.For these model comparisons, no below-sea-level overdeepenings develop.While this is a simplification of the complex feedbacks among the glacial flow, subglacial and englacial hydrology, and the subglacial sediment (Hooke, 1991;Alley et al., 2003), glaciated valleys are produced that fit the expectations that overdeepened glaciers tend to reach critical angles.For the model comparisons, all processes except for glacier sliding and deformation are calculated within the ICE-Cascade model not within the HO nested model.Certain processes, such as avalanching and hillslope processes can, occur in different locations between the HO and SIA model if the ice extent is significantly different.
Glacial erosion is performed by two major mechanisms: abrasion and quarrying (Hallet, 1979(Hallet, , 1996;;Iverson, 1991).These processes operate on spatial and temporal scales that are orders of magnitude shorter than those of the glacial valley, the climate cycle, and the orogen.The erosional processes are often simply upscaled in landscape evolution models (Tomkin, 2003;Herman and Braun, 2008).Following the methods of many existing studies of glacial erosion on orogenic timescales, we use a simple relationship between erosion and sliding.For both models, glacial erosion, ∂z b ∂t , is proportional to the sliding velocity, which comes from modeling and empirical studies of glacial erosion (Harbor et al., 1988;Humphrey and Raymond, 1994), such that where K is a constant that characterizes the erodibility of the subglacial material (Laitakari et al., 1985;MacGregor et al., 2009;Duhnforth et al., 2010) and l is another constant generally equal to 1 (Table 1).While a few recent studies have used more sophisticated rules for erosion (MacGregor et al., 2009;Iverson, 2012), this is the same rule as used in previous ICE-Cascade and other glacially influenced, landscape evolution models (e.g., Braun et al., 1999;Herman and Braun, www.earth-surf-dynam.net/3/153/2015/Earth Surf.Dynam., 3, 153-170, 2015 2008; Kessler et al., 2008;Egholm et al., 2012b;Yanites and Ehlers, 2012).

Shallow ice approximation
The shallow ice approximation (SIA) simplifies full-stress glacial flow by assuming that the ice is significantly wider than it is thick, and the surface slopes are not large, i.e., it is shallow (Fowler and Larson, 1978;Hutter, 1983).In this assumption, all stresses but the simple shear stresses in the direction of ice flow are assumed to be negligible; longitudinal and lateral stresses, including drag, are assumed insignificant.This considerably simplifies how ice flow can be modeled, which is particularly useful when used on orogenic length and timescales.However, the assumptions based on glacial geometry and surface slope, while originally derived for use on large, shallow ice sheets, are not necessarily true when used for alpine glaciers with considerable sliding and where narrow valleys channelize flow and large surface and basal slopes are present (Hutter, 1983;Egholm et al., 2011).Ice thickness, H , is computed from mass-balance equation where M is the mass balance and F is the vertically integrated ice flux.
The ice velocity, u, is simplified to just two dimensions where vertical velocities are deemed negligible, so that u = u î + v ĵ .Each directional component is the sum of two components of glacial motion, which are designated using the SIA in this definition: deformation u sia d and sliding u sia sl for the velocity in the x direction (v is defined in a parallel fashion for the velocity in the y direction).
The deformation component is defined where A, which is pressure-and temperature-dependent but treated as a constant for most of this analysis and is furthered discussed in Sect.4.4, and n come from the Glen flow equation (Cuffey and Paterson, 2010), which relates the strain rate, ε, to the stress, σ , via β is a constriction factor that is dependent upon the gradient of the subglacial topography.In a standard SIA model, β is unity.For this version of ICE-Cascade, the SIA model is modified to incorporate this factor to partially account for the stresses in steep alpine glaciated valleys where the SIA is too simple (Herman and Braun, 2008): The constant, k c , has been evaluated over a 1 km wide glacial valley by Egholm et al. (2011).A value of 1000 was found to have reasonable agreement with a higher-order approximation, with comparisons of sliding and deformation velocities.
Using the constriction factor, overpredictions of erosion rate from the standard SIA were diminished to underpredictions in comparison to the erosion rates from the higher-order approximation (Egholm et al., 2011(Egholm et al., , 2012b)).
The sliding velocity u sia sl is defined in a similar form to the deformation, dependent upon the shear stress at the bed of the glacier, incorporating the sliding flow factor A sl , the constriction factor discussed above β (Table 1), and simple subglacial hydrology as the effective pressure, the difference between the ice overburden pressure and the water pressure, N −P .Water pressure and its change over time are influences on the sliding velocity and has been used in a variety of other glacial erosion and landscape evolution models (e.g., MacGregor et al., 2000;Tomkin, 2003;Herman and Braun, 2008;Kessler et al., 2008;Egholm et al., 2011).However, in order to focus on just the physics of the ice flow, we do not consider variable subglacial hydrology in this model and treat N − P as 80 % of the ice overburden pressure; Sect.4.4 contains more discussion of this choice.

Higher-order physics
While a variety of various higher-order simplifications have been used to represent flow in other models (Pattyn, 1996;Egholm et al., 2011), we opt for a full-stress solution nested into the larger ICE-Cascade framework.Nesting of this modeling within the SIA model is required due to computational considerations, and our analysis is mostly focused on differences between the two model setups over a limited region.Figure 1 shows an example of topography during a glacial maximum with the nested region highlighted (from 100 to 150 km in the x direction and 100 to 160 km in the y direction).In this region, we use the full-stress equations and treat ice as a nonlinear viscous flow.The sliding velocity is calculated for this region and then passed back to ICE-Cascade.
To ensure that the pattern of sliding velocity is smooth, a linear interpolation is used between the sliding velocity at the edges of the nested domain and those in the nested domain.This small region surrounding the nested domain has a width of less than 1 km.
Conservation of mass is given by where u is the 3-D velocity vector u = u î + v ĵ + w k, and is also subdivided into deformational and sliding portions: u = u ho d + u ho sl .Ice flows as an incompressible laminar material.The Glen flow equation, when written in viscosity form for a single stress component, is and the viscosity η is given as where ε is the effective strain rate, the second invariant of the strain rate.A is a flow constant dependent upon the temperature of the ice but treated as a constant in this analysis.
The full strain rate tensor is defined as and the second invariant that is used to define the effective strain rate in the viscosity term, Eq. ( 10), is The shear stress (τ ij ) can be found from Eq. ( 9).With Eqs. ( 9), (10), and (12) combined, the full shear stress acting in the x direction, τ xz , is then defined and this is the term (along with the shear stress in the y direction τ yz ), when defined at the bed of the glacier τ b , upon which the basal sliding is dependent.The sliding velocity is defined in both the x and y directions, similar to Eq. ( 7).For the x direction, the sliding velocity is where A sl is the same constant as in Eq. ( 7), and m is 3 in this case (Table 1).As in the SIA model, the effective pressure (N − P ) is 80 % the ice overburden pressure.

COMSOL Multiphysics
COMSOL Multiphysics has been used for modeling ice flow in 2-D flow-band and flow-line profiles (Johnson and Staiger, 2007;Campbell, 2009;Headley et al., 2012).Here, we modify the steady-state viscous flow module using the Glen flow law.The flow is represented in equation form by the 3-D incompressible laminar flow, and we set the viscosity to be dependent upon the strain rates (Eq.9).The geometry is defined from the bed topography and the ice thickness within the nested domain.The boundary conditions on the edges of the domain are open boundaries.These boundaries are only used within COMSOL simulations, as the velocity distribution over depth is not output back into ICE-Cascade.The top surface is a free surface, and the shear stress at the bottom boundary, along with the basal temperature, is used to calculate the sliding velocity per Eq. ( 14).Within COMSOL, proper meshing is important for solution convergence.In this case, we set the mesh size along the top and bottom ice surface as COMSOL setting "Coarse".From the nonlinear ice flow law, Eq. ( 5), velocity gradients are larger closer to the bed.For more efficient computations, we use nine boundary layers perpendicular to the bed and surface, with decreasing vertical dimension approaching the bed (z direction), using order 10 4 tetrahedral elements; this corresponds to a resolution at the bed spaced between 500 and 800 m in the x and y directions.Figure 2  over a small region (< 1 km wide) bordering the nested region to ensure that the final results are numerically smooth.
The velocities and basal shear stress are interpolated to the node spacing used in the ICE-Cascade model (Herman and Braun, 2008), which is variable to down to 100 m but initially calculated with a spacing of about 1 km.To investigate whether any significant numerical differences had arisen when the models were combined, a preliminary comparison was performed where no erosion took place but the sliding and ice deformation only influenced the glaciated area; no statistical significance in landscape features was found in this comparison, though glaciated area was still varied on a similar scale between the two models.

Results
The ICE-Cascade model outputs used in this study include the glacial erosion rate, subglacial topography, glacier thickness, and sediment thickness; these primary outputs are analyzed for the rest of the study.Figure 1 shows example output of the topography and ice coverage at a glacial maximum.
The size of the nested model region (black box in Fig. 1) was set due to computational limitations in the Stokes flow simulations.The two profiles highlighted in Fig. 1 are used for comparison of the models: one profile is orogen parallel (A-A'), and the other follows a valley profile (B-B').Three separate simulations are performed and compared.The three climate and thermal settings (Simulations 1-3) are further discussed in Sect.2.2, summarized in Table 2, and shown in Fig. 3.For these three simulations, each was performed using both the ICE-Cascade with SIA glacial-flow physics and with the HO nested subdomain.Erosion rate is generally used instead of sliding velocity or basal shear stress because erosion rate not only reflects the patterns of both of these and freezing at the bed of the glacier (Eqs. 1, 7, and 14) but can also be compared to other destructive and constructive processes incorporated into the orogenic model.

Influence of glacial-flow physics on glacial area and volume
Large-scale properties of glaciations can be used to compare the different climate scenarios and the different glacialflow models.In this case, we use the glacier-covered area and the maximum ice thickness for each time step.As expected, colder climates lead to significantly larger glacialcovered areas and persistent thick glaciers (Fig. 4), regardless of glacial-flow model.The glaciations are progressively smaller for each subsequent glaciation even though the temperature amplitude is not changing.
The overall pattern of ice cap growth is similar between the SIA and HO models, and the differences between the higherorder glacial-flow physics and the shallow ice approximation models are generally difficult to see when comparing the full magnitudes of ice-covered area and maximum ice thickness (Fig. 4).However, comparing the percent difference between the two ( A SIA − A HO A SIA × 100 and , where A is the ice-covered area and H * is the maximum thickness, and the denominators are those values averaged over the glaciation center around the 100 kyr) shows significant variation (Fig. 4), reaching a difference of over 30 % for some of these simulations.The percent differences between the two glacialflow models show that even a higher-order nested model that does not cover the full ice-covered domain can have an effect on the fully glaciated area.
At the start of the simulations, the percent differences in the area are small for all the simulations (Fig. 4), and the maximum ice thickness follows a similar pattern (Fig. 5).Ice thickness tends to grow rapidly, reaches a maximum when the maximum area and minimum temperature are reached, and then tapers slightly before decreasing rapidly.As time progresses, changes are exacerbated, and the difference between SIA and HO is more readily apparent.In general, the simulations that contain warmer temperatures and/or have more sliding (i.e., Simulation 3) have a much larger differences in ice thickness and glaciated area due to the HO glacial-flow physics, while Simulation 2, the coldest, has minimal differences.
A null simulation was also performed where the SIA and HO models were used but glacial erosion was turned off.In this case, for both the maximum ice thickness and the glaciated area, differences between the two models were on the order of 5 % or smaller during every glaciation.The differences did not increase as time progressed.This result is presented to emphasize that the differences between the two models for the three simulations come mainly from the models' influence on the glacial erosion rate and are not inherent from the model nesting, model choice, or other erosional processes.

Simulation 1: the influence of glacial-flow physics on glacial erosion pattern
Results and comparisons of the pattern of erosion rate from Simulation 1 are shown in detail in Figs. 5 and 6.Because Simulation 1 uses climate parameters (Fig. 3 and Table 2) that are generally in the middle of the range covered by the simulations, the results from this show the effects of the different physics models without extreme behavior, such as sliding everywhere in Simulation 3 or the extensive frozen bed in Simulation 2. However, since the erosion rate is dependent upon the basal temperature such that, in very cold excursions and at high altitudes, the glacier can be frozen to the bed and those regions are protected from erosion.We look in detail at the pattern of erosion over the area of the nested model (Fig. 1) and compare these patterns averaged over a single glaciation and over the full simulation (400 kyr).
Figure 5 shows the erosion rate averaged over the first full glaciation centered around 100 kyr (grey-shaded region in Figs. 3 and 4) for the nested subdomain, and the corresponding averaged erosion rates over the full simulation are shown in Fig. 6.The patterns of erosion are similar over the two timescales for each of the SIA (Figs. 5a and 6a) and HO (Figs. 5b and 6b) simulations.Within a given model, the erosion rates are more than half as low on the long term, due to the averaging including the interglacial periods when no glaciers are present.On the long term (Fig. 6), the area actively eroded by the glacier is more extensive than over the single glaciation (Fig. 5), as the regions under frozen ice over a single glaciation are not necessarily always under frozen ice over the entire 400 kyr.
Figures 5c and 6c show the differences between the SIA modeled erosion rates and those from the HO model.When comparing the SIA model with the HO model, we focus on the full 400 kyr average (Fig. 6c), as results for the single glaciation average (Fig. 5) follow a similar pattern.In Fig. 6b, the erosion rate in the HO model peaks at over 4 times that of the SIA, yet the bed is frozen over more of the HO model domain.Differences are noticeably larger in regions where the basal temperature significantly varies between the two simulations, i.e., where one model or the other has no sliding (thus no erosion) occurring.While for both glacial physics models there is generally a region of highest erosion rate in the SE corner of the model, the SIA modeled pattern has only one broad region of high erosion rates around x = 150 km and y = 112.The HO modeled erosion rate pattern shows three specific smaller areas of higher erosion rate in the region between x = 140 and 150 km and between y = 110 and 130 km.

Influence of glacial-flow physics on glacial erosion rate in different climates
In order to compare the erosion rates over both the 400 kyr and the single glaciation timescales for all climate simulations, we average the erosion rates over the profiles within the nested domain (Fig. 1). Figure 7 shows the erosion rates for the orogen-parallel (A-A') swath, and Fig. 8 along the valley profile (B-B').Similar to the 2-D contour plots for Simulation 1, the SIA modeled erosion rates (Figs.7a and 8a) are generally larger and more variable than those from the HO model (Figs.7b and 8b).The differences between the two models (Figs.7c and 8c) can be up to almost an order of magnitude greater than the SIA erosion rates, and the differences between the two glacial-flow models are higher for the warmer and wetter simulations.When comparing the differences, those purely related to the climate are also striking.The patterns of erosion within a given glacial-flow model are generally similar in the simulations, and only the magnitudes differ.The differences are largest in the valleys, where ice is thickest and mov- largest erosion rate differences occurring around x = 100 and 132 km (Fig. 7c).

Influence of glacial-flow physics on subglacial topography and sedimentation
Subglacial topography (Figs.9a-b and 10a-b) is composed of not only bedrock but also sediment deposited by proglacial fluvial processes.The effects of the choice of glacial-flow physics can be seen in comparing the topography (bedrock elevation and sediment thickness combined), the bedrock topography, and the sediment layer thickness.These are all related, but there are many differences among the different physical models and the climate.The topography shown in Figs.9a-b and 10a-b is influenced by glacial, fluvial, and hillslope erosion along with sediment that is accumulated when the fluvial system lacks the carrying capacity to transport it.Hillslopes and steep areas are regions of net erosion, particularly seen in the orogen parallel swath (Fig. 9, around x = 117 km and x = 130-140 km), even when glacial erosion rates might be small or nonexistent (Figs. 7 and 8) due to a frozen bed.When comparing the bedrock elevation in the SIA to the HO models (Figs.9d and 10d), the bedrock differences (SIA-HO) have a similar pattern to those of the total topography, though the magnitude is slightly subdued.For the orogen parallel swath (Fig. 9, around x = 120-125 km), significant valley fill only occurs down glacier of the swath, particularly for Simulation 3, where material is eroded everywhere on the glacier, including the regions frozen to the bed in the other simulations.
The differences between the SIA and HO models for total topography, bedrock elevation, and sediment layer thickness (Figs.9c-e and 10c-e) show similar patterns to the erosion rate (Figs.7 and 8).In Figs. 9 and 10, the simulations generally show very similar patterns in both the elevation and the differences between the physical models.Warmer and wetter runs are associated with larger differences between the physical models.Simulation 3 shows the most extreme changes to the topography as well as the most extreme sediment accumulation (Figs.9a-b and 10a-b) and differences between the physical models (Figs.9c-e and 10c-e), whereas Simulation 2 shows the least change except for a large amount of deposition in the SIA model around x = 110 km (Fig. 9e).

Discussion
If a more complex model and a simpler model can be in agreement when the assumptions and approximations in the simpler model are valid, then the simpler model with fewer free parameters is preferred.Model choice also depends upon how this similarity is defined and what the area of interest is, as well as what features or processes are being modeled and what time and length scales are of interest.When looking at a full orogen, it seems that the modified SIA can reproduce features similar to those found in actual orogens.However, in mountainous topography, particularly at sub-polar latitudes, glacier dynamics are influenced by the physical constraints of valleys and fjords and also are a strong control on the erosion and sliding rates.
The different glacial-flow models have an effect on the glacier and the topographic evolution of the orogen, although the magnitude of this effect is variable.Additionally, these effects are dependent upon the climate.When the glacier is mostly frozen (Simulation 2), the physics chosen have only a small influence on the glacial extent and thickness (Fig. 4) and on the topography (Figs. 9 and 10).However, if the glacier has large warm periods (Simulation 1) or is forced to be wet-based (Simulation 3), even if cold temperatures are reached over large periods, then the model choice is considerably more important.
Sliding velocity, which is dependent upon basal shear stress and the basal temperature, is ultimately one of the most important functions in the simulations.As such, erosion rate is used in many of the comparisons and can be considered a stand-in for sliding velocity since our chosen law for erosion rate scales directly with sliding velocity (Eq.1).Comparisons between Simulations 1 and 3, which have the same climate parameters (Fig. 3), emphasize how important the choice of sliding law and reliance on basal temperature are, regardless of glacial-flow model.For example, Figs. 7 and 8 show maximum erosion rates in Simulation 3 of more than double those of Simulation 1. Simulation 1 generally has more than double the erosion rates compared to those of Simulation 2 (Figs.7 and 8), which stresses how important the temperature can be even without extreme erosion laws like in Simulation 3.

Glacial properties and uncertainty in differing climates and over different timescales
Regions like the Gamburtsev Subglacial Mountains underneath the East Antarctic Ice Sheet (Young et al., 2011) illustrate that, no matter how large and thick ice coverage might be, as in Simulation 2 (Fig. 4), if the basal temperature is regularly below freezing, there will be little modification of subglacial topography because there is no sliding.It follows that if the ice is generally below freezing, the glacial-flow model does not tremendously matter if interest is on the evolution of the landscape; in that case, an SIA model could be a reasonable and computationally efficient choice.However, if interest is in ice extent or coverage, which would be important for climate or sea-level-rise studies, then the glacialflow model can be important even for completely cold-based ice when timescales are long enough (Fig. 4, particularly Simulation 2).While it might be predicted that, over multiple glacial-interglacial cycles on geologic timescales, differences between the two glacial-flow physics might be averaged out due to the other erosional processes reshaping the landscape during interglacial periods, that does not appear to .Erosion rates over the A-A' orogen-parallel swath profile (Fig. 1).(a) SIA erosion rates averaged over the 100 kyr glaciation (dashed line; grey-shaded area in Figs. 3 and 4) and over the full simulation (solid line).(b) HO erosion rates averaged over the 100 kyr glaciation (dashed line) and over the full simulation (solid line).(c) Difference between the time-averaged erosion rates (SIA-HO).be the case.Cumulative effects of sediment deposition and transport due to the non-glacial processes in our simulations do not noticeably moderate the final topography when comparing between the glacial-flow physics.
Comparing the simulations among the different climate scenarios and not just between the two physical models allows us to consider the influence of glacial-flow physics versus climate.The climate simulations (Fig. 3) are considerably different, and their effects on the erosional pattern and topographic evolution are substantial: the ice-covered area varies by over a factor of 3 (Fig. 4a), and the erosion rates can vary by over a factor of 2 (Figs.7a-b and 8a-b).Particularly for the orogen-wide properties (ice-covered area and maximum thickness), the variations from the climate are substantially larger than those from the choice of glacial-flow physics model (Fig. 4).However, when the erosion rates and topographic evolution are compared over the swath profiles (Figs.7-10), the differences from the choice of glacial-flow physics model are generally of the same magnitude or larger than those from different climates.These results emphasize that the choice of glacial-flow physics model is less important than the climate if interest is only on larger, orogen-wide properties.For properties like bed topography and erosion rate on the valley scale, the choice of glacial-flow physics can make a more significant difference than even very different climate models.

Evolution of subglacial topography and sediment thickness
The subglacial topography and deposited sediment are important metrics to consider, as in real landscapes these are the relics from previous glaciations, formed through erosion and other geomorphological processes.We evaluate not just how the glacial-flow physics model influences the erosion of topography but also how this topography and the eroded sediment work within the other geomorphological systems.Fluvial erosion is an important mechanism, and the fluvial network is also responsible for the transport and redistribution of glacially eroded sediment once it has left the toe of the glacier (Hallet et al., 1996;Alley et al., 2003).In some cases, the existing rivers do not have the carrying capacity to support the evacuation of all the available eroded material.This material impacts the landscape not only by protecting the bedrock from further erosion but also by impacting the glacial flow, subglacial hydrology, and erosion patterns (Humphrey and Raymond, 1994;Hallet et al., 1996).The sediment thicknesses produced vary slightly between the SIA and the HO models.Along the swath profiles, the pattern of deposition varies between the physical models.Generally, the HO model leads to no change or smaller sediment thicknesses by 5-20 m (Figs.9c and 10c).1).(a) SIA erosion rates averaged over the 100 kyr glaciation (dashed line; grey-shaded areas in Figs. 3 and 4) and over the full simulation (solid line).(b) HO erosion rates averaged over the 100 kyr glaciation (dashed line) and over the full simulation (solid line).(c) Difference between the time-averaged erosion rates (SIA-HO).
in the orogen parallel swath and significantly less (almost 20 m) deposition around x = 120-125 km in the same swath.Along the valley profile B-B' (Fig. 10) there is also considerably less deposition for Simulation 3 in the HO model lower in the profile (Fig. 10c), despite significantly more erosion occurring upstream.These variations in amount of sedimentary fill have implications for the structure of future drainage areas, glacial flow, and isostasy.

Comparison to other models
This study compliments a variety of other research on the effects of glaciers on developing topography.In general, this study compares well with models that also test different ice physics models.However, the magnitude of the effects of climate and ice physics can also be compared to the role of other processes or properties, glacial or otherwise, in landscape evolution models, such as choice of subglacial hydrology, erosion law, and initial conditions.
As mentioned briefly in the model setup section, the sliding rule used is a simple one that does not incorporate subglacial fluvial water pressures or pressure changes.For both sliding and erosion, subglacial water has increasingly been shown to be important (Clarke, 2005;Cohen et al., 2006;Bartholomaus et al., 2008).While much current research, modeling or field-based, is focused on understanding the dy-namics of the subglacial fluvial system on the short-term scale, how this system can be meaningfully scaled up to glaciological or geological timescales is still an open question.In previous work, hydrology has been shown to make a marked difference in the development of topography.Systems with a coupled, dynamic hydrologic system in general led to more sliding and therefore more erosion (Egholm et al., 2012a).Feedbacks developed among the water, topography, and glacier, and more erosion occurred over a shorter period of time than in a control study.In addition, hydrology was found to be linked to the development of topography associated with previously glaciated regions, including hanging valleys and overdeepenings (Egholm et al., 2012a).The effects of the hydrologic system were in general on the order of or larger than most effects discussed within this current study.To isolate just the influence of ice physic, the hydrology has been kept simple.
The choice of erosion law, Eq. ( 1), is based upon another assumption, and various other erosion laws tie the erosion rate to other powers of the sliding velocity (Hallet, 1979;Iverson, 1995;MacGregor et al., 2009), the ice flux (Kessler et al., 2008), or the basal shear stress (Pollard and Deconto, 2007), and a different choice would influence the patterns of erosion and the locations of maximum erosion rate.Comparing the magnitude of effects of erosion laws on developing topography can be difficult for a number of reasons: scaling factors and constants might vary; the ice physics, which can add further complications, determine the basal shear stress; and subglacial hydrology is still quite complicated.Mac-Gregor et al. (2009) determined that the use of a composite erosion law, incorporating both quarrying and abrasion, in a varying climate led to the creation of more rugged terrain with higher-frequency roughness and a deeper cirque (by over 40 %) due to the different erosion laws focusing erosion in different regions along the glacier valley length compared to using the empirical law (Eq.1).In terms of magnitude of change, the use of different erosion laws and different dependencies on hydrology within a non-evolving flow-line model led to differences in erosion magnitude of over 20 % along with shifting of erosion peaks by tens of kilometers (Headley et al., 2012).Over long timescales, this could significantly influence the development of the valley.This variation in erosion rate due to the incorporation of different or composite erosion laws could certainty lead to variations in developing topography larger than those just from changing the ice physics within a fixed glacial erosion law.
Recent work has also emphasized the importance of initial conditions on landscape evolution in glaciated regions (Pedersen et al., 2013).By using topography with more or less relief or that which has been preconditioned by rivers versus glaciers, the final landscape can significantly vary.The initial topography within this study is an equilibrium land- scape generated using only fluvial and hillslope processes (Yanites and Ehlers, 2012).Using glacial steady-state landscapes or non-equilibrium landscapes might lead to different landforms emerging due to the changing focus of glacier erosion.Glaciers occupying fluvial landscapes tend to erode quickly at high altitudes, creating basins within which future glaciers form but which also lower over time so as to create a negative feedback among the glacial erosion and extent (Pedersen et al., 2013).This negative feedback is likely the same as seen in Fig. 4; as the landscape developments, the glaciations become smaller and smaller.This effect might have been further exacerbated had glacial topography been the initial condition of the model comparisons.However, this negative feedback exists within this current study regardless of model choice, so the influences of ice physics would likely still stand out on their own.

Model caveats and limitations
The SIA and HO models used in this study have several caveats and limitations worth mention.First, this study provides only minimal estimates for the effects of the higherorder flow model on the topographic and glacial evolution, particularly over the full domain.Since the higher-order model is nested within the standard model, not all regions of the glacier and the bed feel the effects of the flow, though the nested region covers the ice divide and spans multiple valleys (Fig. 1b and c).However, we suspect that incorporating all of the glacial ice would lead to larger differences in how topography develops.The full orogen comparisons (Fig. 4) might then be viewed as conservative estimates as to how much the HO model can influence the glacier extent and thickness.Many of the simplifications that are incorporated herein act to further emphasize the singular role of choice of glacial physics without incorporating other complicated feedbacks.
The topography falls within the limits considered generally outside of the scope of the SIA in other studies (Egholm et al., 2011).For topography used herein, there is on average about 300 m of valley-ridge relief in the mountain range, with a maximum elevation around 1600 m.The smaller valleys where glaciers flow are generally between 5 and 7 km wide.Ice thickness reaches a maximum of over 1 km and during glaciations is generally hundreds of meters thick.This range of thickness is on the same scale as, or only an order of magnitude smaller than, the valley width.One way of determining whether the SIA is applicable can be through comparison of aspect ratios, such as ice thickness to width or ice thickness to characteristic length (Hutter, 1983).The ratios, particularly of ice thickness to valley width, in this topography generally fall outside the range of usage of the SIA, where the aspect ratio should be 0.1 or lower.In addition, the topography is sufficiently steep beyond the reaches of the SIA.The topography used in this analysis can overall be described as a triangular, east-west-trending mountain range with a large shelf to the south and a smaller shelf to the north, with an average slope around 4 • on the southern slope and a steeper 5 • for the northern slope.The slope at the front of the range, representing the continental shelf, is set at around 0.05 • (a gradient of 0.001).On the valley scale and smaller, the topography has a maximum slope around 20 • , with over 50 % of the topography steeper than 1 • .As the model evolves, the orogen-scale relief and slopes do not vary significantly on the timescales under consideration.
One simplification in this comparison is that of the choice of non-thermomechanically coupled ice flow.The constant A from the Glen flow Law, Eq. ( 5), is treated as constant for all simulations, though this term is actually temperature and pressure dependent.Preliminary comparisons between thermomechanically coupled and uniform temperature simulations found long-term differences between the simulations to be minimal, an order of magnitude less than the differences between the simulations with different glacial physics.
Another erosional feedback that would influence the evolution of the glacier and the topography would come out of the choice of erodibility, K in Eq. ( 1).In this model, all bedrock material is treated with the same erodibility, with no accounting for fracture mechanics or different rock types.In addition, sediment that is eroded then deposited is reincorporated into the model with the same erodibility as the original bedrock.A variety of research related to modeling and field data of different erosional systems has shown that the rock erodibility, whether this varies due to fracturing, rheology, or composition, can have a significant effect on the ultimate form of the landscape (Duhnforth et al., 2010;Ward et al., 2012).
On geologic timescales, glaciations can significantly change not just the extent of a future glaciation but also how future glaciers grow and retreat, which further influences the future patterns of erosion (Pedersen and Egholm, 2013).Similarly, in this study, the differences between the physical models become more evident the longer the topography develops, as multiple glaciations erode the landscape and the fluvial and hillslope processes do not completely overprint these glacial signals.As models move to higher-order physical representations, the question of how well a given model represents reality needs to be further addressed.While this current study focuses on model comparisons, examples of regional hypsometry, patterns of glacial erosion and ice thickness, and patterns of deposition might be used to better determine whether more sophisticated models better reproduce more realistic topography.

Conclusions
This research shows a comparison of a nested HO glacialflow model with a SIA glacial-flow model within the ICE-Cascade landscape evolution model.The simulations incorporate constant tectonic uplift rate, orographic precipitation, hillslope erosion, fluvial erosion, and sedimentation in addition to the glacial erosion processes.Using multiple climate simulations, the effect of glacial-flow physics is evaluated over 400 kyr of topographic evolution.In general, the glacial-flow model choice makes a difference in the development of a glaciated landscape over the long-term, which corresponds to what other studies have shown in simulations without fluctuating temperatures (Egholm et al., 2012b).We have evaluated a variety of properties of landscape evolution in a glaciated orogen, from the ice-covered area and ice thickness to the bed topography and sediment thickness.
Two major conclusions can be reached.First, as expected, the climate, and particularly whether the glacier is mostly cold-based or mostly sliding, has a large influence on glacial development; a few degrees' difference in the minimum sea level can more than triple the glaciated area over an orogen.Changing the climate parameters can lead to large variations in average erosion rates and topography after multiple glaciations have occurred, generally by factors of 3 to 4 between warm and wet glaciers and cold-based glaciers .Second, though these climate influences are large, comparing a HO glacial-flow model to a SIA model show that the choice of model can have as large, if not larger, an influence on the developing orogeny as climate alone for glaciers with more warm beds and less influence for colder glaciers.As topography develops, even with other processes (fluvial, hillslope) dominating over large time periods, the deviation between SIA and HO grows over time to make a difference of over 30 % for completely wet-based glaciers (Fig. 4).The sub-glacial topography is similarly affected by the incorporation of the HO model, as long-term erosion rates can vary by over an order of magnitude between the HO and the SIA models , and the subglacial topography can vary by over 100 m between the physics models over 400 kyr .Differences between the models for the colder climate are significantly subdued.Therefore, for modeling studies of landscape dynamics in cold (generally below freezing) climates, glacier physics based upon the SIA approximation would be the most efficient choice; any sacrifice of more realistic glacier dynamics is not significant.However, over orogenic timescales and relatively warm climates, or if interests are on glacial extent and ice thickness regardless of climate, the choice of higher-order glacial-flow models should be considered, as, compared to the SIA model, such models can lead to large differences in the glacial coverage, developing topography, sediment deposition, and averaged erosion rates.

Figure 1 .
Figure 1.(a) Model domain with ice coverage from glacial maximum for Simulation 2 with hillshade topography at T = 100 kyr.Glacial maxima for the other simulations generally follow the same form.The black box outlines region of nested, higher-order physics domain, used in Figs. 5 and 6.The shaded region (A-A') shows the location of the orogen-parallel swath used for comparisons in Figs.7 and 9. Line B-B' gives the valley profile used for comparisons in Figs. 8 and 10.(b) Bedrock topography and ice surface at glacial maxima of T = 100 kyr and initial bedrock topography along A-A'.(c) Bedrock topography and ice surface at the same glacial maxima and initial bedrock topography along B-B'.
Figure 2. Graphical representation of the mesh used within the HO nested model.The width and depth correspond to those of the black box in Fig. 1.The height corresponds to the ice thickness, and the bottom edge is at the elevation of the topography.The mesh is explained further in Sect.2.5.
shows an example of the mesh used in COMSOL for the nested model.When the sliding and surface velocities and basal shear stress are input back into the ICE-Cascade model, the velocity profile of both the HO region and the surrounding SIA region is linearly interpolated www.earth-surf-dynam.net/3/153/2015/Earth Surf.Dynam., 3, 153-170, 2015

Figure 3 .
Figure 3. Climate variations over time for the different simulations.Over three complete glaciations are simulated.Simulation 3 has the same climate as Simulation 1 and is not specifically shown.The grey box shows the single glaciation used in the analysis presented in Figs. 5, 7, and 8.

Figure 4 .
Figure 4. Variations in glaciated area and maximum thickness over 400 kyr.The colors correspond to the simulations in Table 2 and described in Sect.2.2.For (a) and (c), solid lines indicate the SIA model, and dashed lines the HO.(a) The area is defined at each time step by the area covered by ice greater than 10 m thick.Solid lines indicate the SIA model, and dashed lines the HO.(b) The percent difference between the SIA and HO.(c) The maximum ice thickness is defined at each time step when ice greater than 10 m thick exists.(d) The percent difference between the SIA and HO.

Figure 5 .
Figure 5. Pattern of erosion rates in nested region (Figs. 1 and 2) averaged over a single glaciation (grey-shaded areas in Figs. 3 and 4).The grey-scale topographic lines show the elevation of the underlying topography.(a) SIA erosion rates, where regions of no erosion are shown in white.(b) HO erosion rates, where regions of no erosion are shown in white.(c) Difference between the timeaveraged erosion rates (SIA-HO).

Figure 6 .
Figure 6.Pattern of erosion rates in nested region (Figs. 1 and 2) averaged over the full simulation.The grey-scale topographic lines show the elevation of the underlying topography.(a) SIA erosion rates, where regions of no erosion are shown in white.(b) HO erosion rates, where regions of no erosion are shown in white.(c) Difference between the time-averaged erosion rates (SIA-HO).
Figure7.Erosion rates over the A-A' orogen-parallel swath profile (Fig.1).(a) SIA erosion rates averaged over the 100 kyr glaciation (dashed line; grey-shaded area in Figs.3 and 4) and over the full simulation (solid line).(b) HO erosion rates averaged over the 100 kyr glaciation (dashed line) and over the full simulation (solid line).(c) Difference between the time-averaged erosion rates (SIA-HO).

Figure 8 .
Figure8.Erosion rates following the B-B' valley profile (Fig.1).(a) SIA erosion rates averaged over the 100 kyr glaciation (dashed line; grey-shaded areas in Figs.3 and 4) and over the full simulation (solid line).(b) HO erosion rates averaged over the 100 kyr glaciation (dashed line) and over the full simulation (solid line).(c) Difference between the time-averaged erosion rates (SIA-HO).

Figure 9 .
Figure 9. (a-c) Topographic swath profiles over the A-A' orogen-parallel swath (Fig. 1).Solid lines indicate the ice-bed topography contact, and the dashed lines indicate the bedrock elevation (sediment is allowed to accumulate in ICE-Cascade).(a) SIA topographic profiles at the end of the full simulation.(b) HO topographic profiles at the end of the full simulation.(c) Differences between the topographic profiles erosion rates (SIA-HO).(d) Differences between the bedrock topography (SIA-HO).(e) Differences between the sediment thicknesses (SIA-HO).

Figure 10 .
Figure 10.(a-c) Topographic profiles following the B-B' valley profile (Fig. 1).Solid lines indicate the ice-bed topography contact, and the dashed lines indicate the bedrock elevation (sediment is allowed to accumulate in ICE-Cascade).(a) SIA topographic profiles at the end of the full simulation.(b) HO topographic profiles at the end of the full simulation.(c) Differences between the topographic profiles erosion rates (SIA-HO).(d) Differences between the bedrock topography (SIA-HO).(e) Differences between the sediment thicknesses (SIA-HO).

Table 1 .
Landscape evolution and orographic precipitation model parameters.

Table 2 .
Landscape evolution and orographic precipitation model parameters.