Basal shear stress under alpine glaciers: insights from experiments using the iSOSIA and Elmer/Ice models

Shear stress at the base of glaciers exerts a significant control on basal sliding and hence also glacial erosion in arctic and high-altitude areas. However, the inaccessible nature of glacial beds complicates empirical studies of basal shear stress, and little is therefore known of its spatial and temporal distribution. In this study we seek to improve our understanding of basal shear stress using a higher-order numerical ice model (iSOSIA). In order to test the validity of the higher-order model, we first compare the detailed distribution of basal shear stress in iSOSIA and in a three-dimensional full-Stokes model (Elmer/Ice). We find that iSOSIA and Elmer/Ice predict similar first-order stress and velocity patterns, and that differences are restricted to local variations at length scales of the order of the grid resolution. In addition, we find that subglacial shear stress is relatively uniform and insensitive to subtle changes in local topographic relief. Following the initial comparison studies, we use iSOSIA to investigate changes in basal shear stress as a result of landscape evolution by glacial erosion. The experiments with landscape evolution show that subglacial shear stress decreases as glacial erosion transforms preglacial V-shaped valleys into U-shaped troughs. These findings support the hypothesis that glacial erosion is most efficient in the early stages of glacial landscape development.


Introduction
The widespread late-Cenozoic glaciations produced distinctive glacial landforms in many mid-to high-latitude mountain ranges (e.g.Penck, 1905;Sugden and John, 1976).The glacial landforms include U-shaped valleys, bowl-shaped cirques, hanging valleys, and truncated spurs.The consistent geometry of these landforms and the associated non-fractal spatial scales show clear links to the dynamics of viscous flow (Evans and McClean, 1995;Pelletier et al., 2010), which indicates that subglacial dynamics must be of first-order importance to landscape evolution (e.g.Harbor et al., 1988;Anderson et al., 2006).However, measures of subglacial dynamics, such as basal shear and normal stress, are inherently difficult to obtain owing to the general inaccessibility of the subglacial environment.
Only a few studies have measured sliding velocity and basal stress directly; examples include the Glacier d'Argentière in the French Alps (Boulton et al., 1979) and under Engabreen in Norway (Cohen et al., 2000(Cohen et al., , 2005;;Iverson et al., 2003).These studies measured regional shear stress values between 0.1 and 0.3 MPa.However, interpretations from these studies are complicated by their limited spatial and temporal extent, and by local heterogeneity such as the presence of cavities that might concentrate stress at much higher values.It is therefore not possible to investigate catchment-wide variations in shear stress from these empirical studies.Knowledge of spatial and temporal variations in subglacial dynamics therefore relies mostly on inversion of geophysical data (e.g.Joughin et al., 2006Joughin et al., , 2012;;Habermann et al., 2013;Morlighem et al., 2013).Despite several complications in such studies (Joughin et al., 2004;Gudmundsson and Raymond, 2008;Habermann et al., 2012), and very different subglacial settings, these studies also find basal shear stress of the order of 0.1-0.4MPa.
Although often hidden by results focussing on subglacial sliding rate, basal shear stress is an important underlying factor for scaling glacial erosion.Erosion rate is commonly assumed to scale with either basal sliding speed (e.g.Oerlemans, 1984;Harbor et al., 1988;Braun et al., 1999;Tomkin, 2009;Egholm et al., 2009;Herman et al., 2011) or ice discharge (e.g.MacGregor et al., 2000MacGregor et al., , 2009;;Anderson et al., 2006;Kessler et al., 2008), and both depend on subglacial stress through sliding relations.Resolving variations in basal stress under glaciers is therefore important for modelling and understanding patterns of glacial erosion.
Ice motion can be computed using the Stokes equations (Stokes, 1845), which balance the stress components in the ice under the assumption of negligible inertia.Solving the full set of Stokes equations is a computationally demanding task, and most applications therefore use computationally efficient shallow ice approximations (Mahaffy, 1976;Hutter, 1983;Blatter, 1995;Baral et al., 2001;Pattyn, 2003;Egholm et al., 2011).However, it is well known that the accuracy of these approximations depends strongly on the aspect ratio of the ice (ice thickness vs. horizontal extent), the bed slope, and horizontal gradients in ice velocity (Hutter, 1983;Baral et al., 2001).
As an end-member approximation, the zeroth-order shallow ice approximation (SIA) is computationally very efficient, but the approximation is only considered valid for the interior parts of large ice sheets where ice surface gradients are small and smoothly varying (Hutter, 1983;Le Meur et al., 2004;Hindmarsh, 2004).The limitation of SIA models arises mainly because the approximation ignores spatial stress gradients that provide regional coupling of ice flow across a glacier.The latter drawback has led to an increased use of higher-order shallow ice models (HOM), which are considered more accurate in cases where ice velocity vary over relatively short distances (e.g.Pattyn, 2003;Hindmarsh, 2004;Egholm et al., 2011).However, the precise relationship between the aspect ratio of the ice and the accuracy of the shallow ice approximations is only vaguely defined.As a rule of thumb, the aspect ratio should be very small (< 10 −2 ) for a zeroth-order approximation like SIA, while it may be higher (up to 1) for a second-order shallow ice approximation (Baral et al., 2001).Thus, although the higher-order ice dynamics of HOMs should increase accuracy compared to SIA models in steep landscapes, they too will be challenged for example when bed slopes increase beyond a certain limit.These limitations and their implications have received little attention in alpine settings, despite being of prime importance to a number of areas in glaciology and landscape evolution.
Existing benchmark studies have compared results from different models (SIA to full-Stokes models) (Hubbard, 2000;Le Meur et al., 2004;Hindmarsh, 2004;Pattyn et al., 2008;Ahlkrona et al., 2013), but all have focussed on simple descriptions of three-dimensional glacial landforms, often formulated by mathematical functions.Using data from the Haig Glacier in the Canadian Rocky Mountains, Adhikari et al. (2013) investigated the effects of higher-order dynamics for the future glacial evolution.Owing to the overdeepened bed, higher-order effects were suppressed as geometric constraints limited the horizontal glacial flow.In a recent study, Headley and Ehlers (2015) compared two glacial models (a SIA model and a three-dimensional full-Stokes model) in a realistic landscape and found marked differences between models.As it is vital for predictions of ice flow and subglacial erosion to resolve subglacial stress accurately, we performed new comparison experiments on a synthetic but realistic three-dimensional landscape using both the iSOSIA higher-order model and the Elmer/Ice full-Stokes model (Sect.2.3).While this setup prevents analytical solution of the Stokes equations, it allows us to compare the iSOSIA approximation to a full-Stokes computational model in a realistic setting under different scales of relief (Sects.3.1 and 3.2).
In subsequent experiments, the same landscape provides the basis for iSOSIA experiments that combine subglacial erosion with different models for basal sliding (Sect.3.3).These final experiments are designed to explore long-term feedbacks between landscape evolution and subglacial dynamics.

Methods
In the following we introduce the ice models used in this study, along with technical details on experimental setup and model comparison.

Elmer/Ice
The Elmer multi-physics software package (www.csc.fi/elmer) provides a finite-element framework for modelling both linear and non-linear three-dimensional flow problems.The Elmer software is developed at CSC in Finland with collaborators around the world, and is published under a GNU Public License (GPL).A special edition of Elmer, named Elmer/Ice, is available with algorithms designed especially for problems related to ice flow (Gagliardini et al., 2013).
Elmer/Ice provides a highly accurate description of glacial dynamics by solving the full set of Stokes equations in three dimensions.However, the high degree of accuracy comes with a high computational demand.Elmer is developed to run very efficiently in parallel (Gagliardini et al., 2013) to reduce computation time, but the computations performed here still required 2 to 3 orders of magnitude more time than the corresponding SIA and iSOSIA simulations.Owing to the high computational demand, we only use Elmer/Ice to perform steady-state simulations without erosion.

iSOSIA
iSOSIA was developed specifically for modelling glacial landscape evolution (Egholm et al., 2011).The ice model includes all stress components of the Stokes equations.However, by using a second-order shallow ice approximation (Baral et al., 2001) iSOSIA represents a computationally efficient alternative to full-Stokes models.The main limiting assumption in iSOSIA is that horizontal, longitudinal, and transverse stress components are not allowed to vary with depth in the ice.This assumption facilitates analytical depth integration of velocities, and iSOSIA is hence a depthintegrated two-dimensional model.
The iSOSIA equations are highly non-linear because components of stress and ice velocity are connected through the non-Newtonian Glen's flow law for ice with a stress exponent of 3. The non-linear equations are relaxed using an iterative red-black finite-difference Gauss-Seidel method (Briggs et al., 2000).iSOSIA was also recently ported to graphical processing units with increased computational efficiency (Braedstrup et al., 2014).

Experimental Setup
The first two experiments are designed to compare stress and velocity components from iSOSIA to those from Elmer/Ice.The objective is to test how well iSOSIA and Elmer/Ice agree on spatial variations in basal shear stress across gradients in topographic relief, ice thickness, and flow rate.The third experiment is used to study how patterns of basal shear stress and sliding evolve when topography change due to subglacial erosion.All three experiments are performed on a synthetic topography generated using a fluvial landscapeevolution model based on stream-power erosion (Fig. 1a; Braun and Sambridge, 1997).This provides a particularly convenient setup where the uppermost drainage divide follows the grid boundaries, avoiding ice flow out of the model domain.The fluvial landscape has V-shaped valleys and concave longitudinal valley profiles that drain the landscape from a maximum elevation of 2500 m above sea level down to 0 m (Fig. 1a).The computational grid is 20 × 40 km, consisting of 100 × 200 cells (i.e.200 m resolution).
Ice thickness is time-integrated using the continuity equation, where H is ice thickness, t is time, q is ice flux, and M is the rate of ice accumulation/ablation.Accumulation and ablation are modelled as a simple linear function of atmospheric temperature: www.earth-surf-dynam.net/4/159/2016/Earth Surf.Dynam., 4, 159-174, 2016 Atmospheric lapse rate 6.0 where is the atmospheric temperature.T sl is the sea-level temperature, dT h is the lapse rate, and h is bedrock elevation above sea level.m acc is the accumulation gradient and m abl is the ablation gradient.All values are listed in Table 1.Experiments 1 and 2 assume steady state, and use the continuity equation only to construct the steady ice-thickness configuration (Fig. 1b).In experiment 3 the continuity equation is used to update ice thickness throughout transient simulations.However, in order to avoid that feedbacks between mass balance, ice thickness, and topography influence the subglacial stress distribution, we use initial bed elevation in the mass-balance function, and we furthermore fix the mass balance in time and ignore the influence of topographical change by erosion on accumulation and ablation.This invariant mass-balance function prevents that secondary effects related to mass balance mask the differences in stress caused by different sliding and erosion laws.
Ice creep and basal sliding contribute to the ice flux vector, q, in Eq. ( 1).The rate of ice creep is governed by Glen's flow law: where εij is the deviatoric strain rate tensor and s ij is the deviatoric stress tensor.A and n are ice flow parameters (Table 1), and τ e is the effective stress: iSOSIA uses a depth-integrated version of Glen's flow law to compute the depth-averaged flow velocities (Egholm et al., 2011).
The iSOSIA-Elmer/Ice benchmarking experiments 1 and 2 use a simple Weertman sliding relation to relate basal shear stress to the rate of subglacial sliding.The Weertman relation was the only sliding model that successfully converged in Elmer/Ice with our setup.In experiment 3, however, we combine iSOSIA with two additional sliding relations to examine the general sensitivity of subglacial stress to first-order assumptions on basal sliding velocity.The three sliding models are all represented by relations between basal shear stress and sliding velocity: (Weertman, 1957), (6) Empirical sliding: τ n s = u s N/C e (Budd et al., 1979), ( 7) (Schoof, 2005;Gagliardini et al., 2007).
Here τ s is basal shear stress; C w , C e , and C c are sliding coefficients specific to each individual relation (Table 1); u s is basal sliding velocity; and λ 0 is a constant defining the overall bed geometry in the Coulomb-friction sliding model (Schoof, 2005;Gagliardini et al., 2007).The effective pressure, N = t n − p w , is the difference between the ice-bed normal stress, t n , and water pressure, p w .The two latter sliding relations (the empirical sliding model and Coulomb friction), which are both used in experiment 3, depend on effective pressure and hence subglacial water pressure.However, in order to focus on first-order correlations between topography and subglacial shear stress, we simplify the influence of hydrology and initially assume that water pressure is everywhere 80 % of the ice overburden pressure, p w = 0.8ρ i gH .In the final experiment we combine two additional flotation fractions of 70 and 90 % with the Coulomb-friction slid-ing relation to test the influence of water pressure.We note, however, that more complex distributions of meltwater pressure may potentially affect patterns of subglacial shear stress through the influence of sliding (e.g.Flowers and Clarke, 2002;Werder et al., 2013;Beaud et al., 2014).However, such effects are beyond the scope of the present study.
In experiment 3 the glacier erodes its bed according to the following sliding-based erosion law: where ė is erosion rate perpendicular to the bed, K a is the erosion constant (Table 1), m is the erosion exponent, and u s is sliding rate.Parameters governing subglacial erosion through abrasion and sliding are still being debated, and it is particularly relevant to question how well the sliding-based law represents subglacial quarrying (Iverson, 2012).However, sliding-based erosion laws have been shown by models to produce realistic glacial landforms (Harbor, 1992;Seddik et al., 2005;Pedersen et al., 2014), which is why we use one here to study how transformation of a landscape from fluvial style to glacial style influences the patterns of subglacial shear stress.Our use of the above erosion law is thus motivated more by phenomenological arguments (i.e. the erosion law leads to realistic glacial landforms) than empirical evidence.We do however perform all erosion experiments with both a linear (m = 1) and non-linear (m = 2) model in order to evaluate the robustness of our conclusions.

Comparing the output of iSOSIA and Elmer/Ice
To ensure comparability between results produced by iSOSIA and Elmer/Ice, both models operate on the same synthetic input topography, represented by a rectangular grid with specified bed elevation in each grid cell.The iSOSIA solver operates directly on this two-dimensional grid, whereas for Elmer/Ice the two-dimensional grid is extruded to a full three-dimensional mesh with five vertical levels spanning the thickness of the ice.This gridding approach ensures that both models use exactly the same topographic input and mesh topology, except for Elmer/Ice having the additional vertical layering.In order to compare ice dynamics on exactly the same ice configuration, a steady-state ice distribution is generated using iSOSIA and subsequently used by both models.A free-slip boundary condition is implemented along grid edges, and isothermal conditions are assumed everywhere in the grid.Both models compute basal shear stress, τ s , from the Cauchy stress tensor, σ b , at the bed: where σ n is the stress vector perpendicular to the bed, and n b is the normal vector at the bed: We also compare the longitudinal and transverse stress components, s xx , s yy , and s xy from both models.However, since these are only computed by iSOSIA in a depthaveraged version, we need to depth-average the horizontal stress also from Elmer/Ice.We obtain the depth-averaging using the following function: and similarly for s yy and s xy .z is depth below the ice surface and H is local ice thickness.We note that for Elmer/Ice the depth-averaged stress components do not enter the Cauchy stress tensor in Eq. ( 10), which is instead constructed from the local stress state at the base of the ice.
In order to provide a frame of reference, we also compare the basal shear stress of iSOSIA and Elmer/Ice to the driving-stress approximation, which is used as a proxy for basal shear stress in zeroth-order shallow ice approximations (Cuffey and Paterson, 2010).The driving stress is computed as where ρ i is ice density, g is gravitational acceleration, H is ice thickness, and h is the elevation of the ice surface.

Experiment 1 -comparing steady-state solutions
The first two experiments are used to benchmark stress and velocity components from iSOSIA against those from Elmer/Ice.The steady-state ice configuration, which is first computed by iSOSIA and then used as input for both models, includes a main trunk glacier fed by several smaller tributary glaciers (Fig. 1b).Ice thickness reaches a maximum of 700 m in the main valley, and thins towards the glacier front and upwards in the tributaries.The depth-averaged creep velocity is highest where the ice is thickest in the main valley, reaching levels of 120 m a −1 (Fig. 1d).Basal sliding speed is high in the main valley and in the steeper parts of the high tributaries (Fig. 1c).
In experiment 1, the spatial distribution of stress is characterised by similar large-scale patterns in iSOSIA and  Elmer/Ice (Fig. 2a).The components of horizontal normal stress, s xx and s yy , are generally positive at high elevations, which reflects an overall extensional stress state in the accumulation zones.In the trunk valley at lower elevations, both stress components are in places negative (compressive) due to local deceleration of the ice.The latter tendency is, however, clearly affected by the details of the bed topography.The horizontal shear stress, s xy , is large, although of opposite sign, along both sides of the main valley due to a strong velocity gradient perpendicular to the main flow direction.
Differences between horizontal differential stress (s xx , s yy , and s xy ) in Elmer/Ice and iSOSIA are in general below ±0.03 MPa in tributary valleys and ±0.01 MPa in the main valley (Fig. 2a, right column).
The basal shear stress is up to 0.2 MPa under the ice in the main trunk valley and near the tributary headwalls.Between these areas, basal shear stress is rather uniform at levels around 0.1 MPa.Differences between Elmer/Ice and iSOSIA are of order 0.025 MPa but up to 0.05 MPa in a few areas (mostly along ice margins).
Sliding velocities are also similar in both models: ∼ 40 m yr −1 in the trunk valley and around 20 m yr −1 in the tributaries.However, the Elmer/Ice solution has areas with high-frequency variations in basal sliding and shear stress, which are absent in the iSOSIA result.These areas have larger differences in sliding velocity between the two models (Fig. 2c, right column).
To aid the comparison between iSOSIA and Elmer/Ice we extract stress and velocity components along two profiles www.earth-surf-dynam.net/4/159/2016/Earth Surf.Dynam., 4, 159-174, 2016 in the transverse and longitudinal directions of the valley (Fig. 3).Profile A-B runs directly across the main valley, while profile C-D starts at high elevation and follows the ice drainage along the main valley down to sea level.The three horizontal stress components (s xx , s yy , and s xy ) are all of the order of ±0.04 MPa along the profiles, but vary in ways that reflect the bed topography.In the transverse direction (profile A-B, left panels of Fig. 3), stress components generally change sign in response to how the velocity components u x and u y vary across the valley (Fig. 1).The basal shear stress along the profile is 2-4 times greater in magnitude than the horizontal stress components, which highlights how basal shear stress dominates the force balance of valley glaciers.
Along the longitudinal profile (c and d, right panels in Fig. 3), the stress components also fluctuate around zero.A clear anti-correlation exists between s xx and s yy , which indicates horizontal pure shear deformation in response to inflow of ice from the tributaries (Fig. 2a).The basal shear stress is remarkably constant along the profile and decreases only slightly up-glacier.This may seem surprising as bed slope increases significantly up-glacier.However, in this case, the effect of bed slope is counteracted by ice thinning.
There are no clear trends in misfit between iSOSIA and Elmer/Ice and the two models generally predict the same patterns and magnitude of stress.Again, the main difference between results is that high-frequency stress variations are slightly larger for Elmer/Ice than for iSOSIA, particularly so for the basal shear stress (Figs.2b and 3).As expected, the SIA driving stress is generally higher, and shows more intense variation, than the basal shear stress for both iSOSIA and Elmer/Ice.

Experiment 2 -the effect of relief
In the second experiment we gradually increase the total relief of the fluvial landscape to test how this influences the consistency between iSOSIA and Elmer/Ice results.Theoretically, increasing the relief should decrease the accuracy of iSOSIA as bed gradients and spatial variations in flow velocity intensify.
We use a simple scaling of the fluvial topography from experiment 1 in order to systematically increase the relief without affecting the drainage patterns (Fig. 4).We then run iSOSIA to a steady-state ice configuration for all amplified topographies and transfer the resulting ice thickness to Elmer/Ice in order to compute stress and velocity components under similar conditions.
When up-scaling relief, the ice-creep velocity increases significantly, the glacier thins, and its front margin advances (Fig. 4a).Because the ice-flow velocity is amplified almost uniformly, the magnitude of the horizontal stress components, which reflect local velocity gradients, also increases in response to the larger relief (Fig. 5).All three stress components still vary around 0 MPa, but the amplitude of the variation increases with relief.The largest response in horizontal stress due to increased relief occurs in the steep highelevation areas near the headwalls.
In contrast to the englacial horizontal stress, basal shear stress is almost unaffected by the increasing relief and the mean value remains around 0.2 MPa for all four situations (Fig. 5).The local deviation from this trend increases a bit from 0.02 to 0.05 MPa as the landscape steepens.
Examining differences between Elmer/Ice and iSOSIA, we note that both models agree on regional stress patterns, and that iSOSIA stress follows the Elmer/Ice solution reasonably well across the range of relief tested here.The regional misfit remains small (< 0.02 MPa) even when maximum relief is 6250 m (Fig. 5c).There are, however, areas where the comparison exposes an increasing misfit (up to 0.1 MPa) between iSOSIA and Elmer/Ice, particularly when focussing on variations at length scales of a few hundred metres.These areas are mainly associated with thin ice and steep ice-surface topography near the glacial terminus or the headwall areas.
Unlike the basal shear stress from iSOSIA and Elmer/Ice, both regional and local variations in SIA driving stress increase significantly with relief (Fig. 5, blue line).While the misfits between Elmer/Ice and iSOSIA are of the order of 0-0.05 MPa for a relief of 5000 m (with spikes up to 0.1 MPa), the misfit between SIA and Elmer/Ice quickly reaches levels well above 0.2 MPa.This misfit is caused by the driving stress' lack of sensitivity to regional velocity variations as well as bed topography.

Experiment 3 -evolution of stress under glacial erosion
After evaluating steady-state solutions of iSOSIA against Elmer/Ice we now investigate the long-term transient evolution of basal shear stress in response to subglacial erosion and landscape development.We only used iSOSIA for this experiment as the computational costs of Elmer/Ice prevent us from running simulations over the thousand-year timescales required for glacial landscape development.The initial topography from experiment 1 was used as input for iSOSIA and was slowly eroded using a sliding-based erosion law (Eq.9).First, we ran the experiment using the Weertman relation for sliding (Eq.6; Weertman, 1957) in combination with a non-linear erosion law (m = 2 in Eq. 9).We find that the V-shaped fluvial valley structure is transformed into a wider and steep-sided U-shaped trough (Fig. 6).This is in agreement with previous studies (Harbor, 1992;Seddik et al., 2005;Egholm et al., 2012).Several other characteristic glacial landforms also appear as a result of glacial erosion, including steep and narrow upper ridges, flattened valley floors, hanging valleys, and truncated spurs (Fig. 6b).
As expected, bed slopes increase in many areas of the landscape, particularly along valley sides and near headwalls (Fig. 7).However, along the longitudinal flowline of   3).The SIA driving stress approximation is also shown in the third row for comparison (blue line).Fourth row shows the ice-surface and bed topography along the same profiles.Note that elevation is indicated on both the left and right axis.
the glacier, bed slopes generally decrease as glacial erosion flattens the valley floor and removes bedrock features that obstruct flow.This development generally causes bed shear stress to decrease in amplitude and become more uniformly distributed under the ice (Fig. 7).The reduction in basal shear stress also decreases sliding velocity as a result of the Weertman sliding relation (Eq.6; Fig. 7).
To test the robustness of this trend we repeated experiment 3 using two additional sliding relations: the empirical relation (Eq.7; Budd et al., 1979) and the Coulomb-friction relation (Eq. 8;Schoof, 2005;Gagliardini et al., 2007).The different sliding coefficients (Table 1) were calibrated to give similar average rates of sliding.The three sliding relations predict slightly different distributions of subglacial shear stress, but all agree on the first-order patterns and magnitudes (Fig. 8).All three relations predict high shear stress in the trunk valley and upper tributaries, and short-scale stress variations that mimic the details of the valley morphology.With increased erosion all sliding relations lead to decreased and more uniformly distributed basal shear stress.This effect is strongest for the Weertman relation but occur for all three relations.As a consequence of the decreasing stress, the spatially averaged sliding velocity also exhibits an overall decrease with erosion in the trunk valley (Fig. 9a).
To further test the robustness of this trend, we repeated the experiments using: (1) a linear erosion law (m = 1 in Eq. 9) in combination with all three sliding laws (Fig. 9b) and (2) the Coulomb-friction sliding model using two alternative flotation fractions (70 and 90 %) to compute effective pressure, N (Fig. 9c).
The similar outcome of these final experiments suggests that decreasing basal shear stress in response to erosion is largely independent of hydrology and the exponent, m, in the erosion equation (Fig. 9).The latter indicates that first-order topographical change is more important than the details of the erosion law.

The benchmarking experiments
In order to estimate the practical utility of iterative higherorder shallow ice approximations, we have compared the results of two different computational methods: iSOSIA and Elmer/Ice.The comparison experiments were designed to reflect a realistic setting of relevance for long-term glacial landscape-evolution studies.However, the realism of the experiments, involving complex topographical variations, also means that we cannot obtain any exact solution for the stress or velocity distributions, and therefore cannot quantify the true accuracy of any of the two computational methods.Instead, the objective of the comparison study is to estimate the difference between stress predicted by the two methods under conditions that are as similar as possible, and we have designed the experiments to meet this criteria.Both methods use the same bed topography and ice distribution, as well as the same horizontal Cartesian grid structure (with an additional vertical grid dimension for Elmer/Ice).We note here that the true accuracy of both methods is expected to depend on grid resolution.For example, the finite-element method in Elmer/Ice allows irregular grid structures that may increase the computational accuracy, e.g.along ice margins where steep ice surface gradients call for increased spatial resolution (Durand et al., 2011).However, a comparison study designed to uncover the difference caused by various approximations to the Stokes equations is only meaningful if similar meshes are used.
The comparison study shows that iSOSIA and Elmer/Ice predict the same overall patterns of stress and velocity.In both models, components of horizontal stress and stress gra- dients interact with flow patterns on a regional scale (i.e.across topographical gradients, which is in strong contrast to the driving stress in SIA models).This highlights, in accordance with previous studies (e.g.Le Meur et al., 2004;Hindmarsh, 2004;Egholm et al., 2011Egholm et al., , 2012)), the benefits of HOM and full-Stokes models over SIA models.The main difference in results from iSOSIA and Elmer/Ice seems to be confined to spatial scales of a few hundred metres (i.e. the grid-cell spacing).In particular, Elmer/Ice includes high-frequency fluctuations in basal shear stress and sliding (Fig. 2), whereas the iSOSIA results appear more smoothly varying.The smoother pattern in iSOSIA is perhaps not surprising when considering the inherent depth integration of horizontal stress.On the other hand, the relatively large differences in sliding velocity between neighbouring grid cells in Elmer/Ice are surprising and cannot readily be explained by variations in bed topography.We therefore ascribe these high-frequency stress variations in Elmer/Ice to the ice-thickness configuration used, which was generated using iSOSIA and therefore in balance with the governing equations of this ice model.Elmer/Ice is not allowed to make adjustments to the ice thickness that would otherwise dampen the local stress variations.
The high-frequency variations in Elmer/Ice are amplified slightly when the total catchment relief increases from 2500 to 7500 m (Fig. 5).On the other hand, the more regional accordance between iSOSIA and Elmer/Ice stress predictions seems almost unaffected by the increasing relief.This is in contrast to the SIA driving stress, which rises with increasing relief.

The evolution of stress in response to erosion
The iSOSIA simulations with erosion (experiment 3) suggest that variations in basal shear stress are generally reduced by the gradual transformation from a fluvial to a glacial topography.This trend can be recognised for all sliding laws tested in this study (Weertman,empirical,Coulomb friction;Fig. 8) and for two different sliding exponents in the erosion law.
The highest basal shear stress is associated with bends in the fluvial channel profile that form interlocking spurs (Fig. 7).These spurs are truncated by glacial erosion and this decreases basal shear stress.In the main valley, glacial erosion thereby efficiently removes obstacles and straightens the path of ice flow.In addition to this, glacial erosion flattens the longitudinal valley profile and widens its cross section, which also contributes to reduced basal shear stress (Harbor, 1992;Seddik et al., 2005).
It is not surprising in the current study that the modelled glacial erosion primarily attacks portions of the glacial bed where basal shear stress is high.Basal shear stress is connected to sliding rate through a sliding relation (Eqs.6-8), which, in turn, is assumed to scale rates of erosion (Eq.9).We note that different rules for glacial erosion, for example the ones based on mechanics of bedrock quarrying (Iverson, 2012), could depend differently on sliding and shear stress.Feedbacks between stress and erosion might be different in such cases.On the other hand, the experiments presented here result in topographic features that resemble well-known glacial landforms, and it seems reasonable that smoother and flatter post-glacial landforms are associated with less drag from the ice.
With increased erosion, the resulting decrease in basal shear stress leads in our experiment to a lowering of sliding rate, and hence slowdown of erosion.This is a direct consequence of the sliding relations used.All of the sliding relations have a power-law scaling between stress and sliding, and sliding must decrease with decreasing shear stress (the power-law scaling of the Coulomb-friction sliding model occurs below an upper limit to the bed's ability to support shear stress).In addition to the sliding rate, basal shear stress associated with the sliding models depends on bed resistance and roughness, which is controlled by parameters C w , C e , C c , and λ 0 in Eqs.(6-8).Here we speculate that if bed resistance decreases more rapidly than the shear stress imposed on the bed by the ice flow, then sliding may possibly accelerate as topography is eroded, in contrast to the results presented here.The bed resistance could decrease if, for example, the bed is smoothed by erosion or if the flatter glacial longitudinal valley profile reduces the meltwater drainage efficiency of the glacier.Our experiment 3 does not show such behaviour, partly because we ignore complex effects of meltwater hydrology, and partly because the sliding parameters representing bed roughness are treated as constants independent of erosion.However, understanding how glacial erosion affects the topographical conditions that promote cavitation on length scales below the current grid resolution of landscape-evolution models may be important for advancing our understanding of feedbacks between glacial dynamics and topographical development.
Because of general lowering of basal shear stress with erosion, our results support the hypothesis that glacial erosion is most efficient in the initial phase of glacial landscape evolution, before landforms are adapted to the new glacial regime (Harbor, 1992;Braun et al., 1999).In general, however, landscape evolution is influenced by several processes not accounted for here, such as a transient climate forcing (Pedersen and Egholm, 2013), changing topography due to a tectonic forcing (e.g.Tomkin and Braun, 2002), periglacial processes acting in concert with glacial erosion (e.g.Egholm et al., 2015), fluvial processes, and mass wasting affecting the landscape, especially during ice-free interglacial periods (e.g.Schlunegger and Hinderer, 2003), and subglacial fracturing in response to high differential stresses (e.g.Leith et al., 2014).

Conclusions
We have investigated and compared the spatial distribution of subglacial shear stress in both a higher-order shallow ice model (iSOSIA) and a full-Stokes model (Elmer/Ice).Using iSOSIA only, we also investigated the temporal evolution of basal shear stress in response to subglacial erosion.In total, we conducted three experiments in order to resolve different aspects of subglacial shear stress.Our findings are as follows: -iSOSIA and Elmer/Ice produce stress and sliding patterns that are largely similar under the conditions tested.
-In the alpine setting used here, basal shear stress seems rather insensitive to increases in overall relief, as reduction in ice thickness counteracts the effects of bed steepening.Thus, increasing total relief by a factor of 3 only produces a small response in basal shear stress.
-Subglacial erosion removes obstacles that give rise to high basal shear stress in the pre-glacial landscape setting.By this, glacial erosion leads to lower and more uniformly distributed basal shear stress.
-Using three different sliding relations and two different erosion laws, we find a stabilising feedback between basal shear stress, sliding, erosion, and topography.This feedback depends, however, on constant sliding coefficients, which in a more realistic setting could be altered by long-term changes to the bed roughness.

Figure 2 .
Figure 2. (a) Comparison of stress components from iSOSIA and Elmer/Ice shown in top view.The first two columns show the depthaveraged stress components s xx , s yy , and s xy .The right-most column shows the difference between iSOSIA and Elmer/Ice results.(b) Basal shear stress for both models.(c) Sliding velocity using a Weertman relation (Eq.6).Ice flow is from right to left.

Figure 3 .
Figure 3. iSOSIA (orange) and Elmer/Ice (green) stress components in MPa along a transverse (left column) and a longitudinal (right column) profile.Upper three rows compare the higher-order depth-averaged horizontal stress components.The fourth row shows the basal shear stress along the same profiles.The SIA driving stress (Eq.15) is also shown for comparison (blue line).The fourth row also shows bed topography (black line) and ice thickness (blue shaded area).Notice that elevation is indicated on the right axis.Bottom left panel shows ice thickness.The positions of the two profiles (A-B and C-D) are shown in the bottom right panel.Note that the bottom panels are map views of the 3-D model (Fig. 1a, b).

Figure 4 .
Figure 4. Bed topography with ice cover (a), creep velocity (b), and basal sliding velocity (c) computed at steady state using iSOSIA for the four scaling factors in experiment 2. The total relief is 3750, 5000, 6250, and 7500 m respectively.

Figure 5 .
Figure 5. Stress components in MPa from iSOSIA (orange) and Elmer/Ice (green) for the increasing topographical relief in experiment 2. Left column of each panel (a-d) shows values along the transverse profile, A-B, while the right column is along the longitudinal profile, C-D (Fig.3).The SIA driving stress approximation is also shown in the third row for comparison (blue line).Fourth row shows the ice-surface and bed topography along the same profiles.Note that elevation is indicated on both the left and right axis.

Figure 6 .Figure 7 .Figure 8 .
Figure 6.(a) Preglacial landscape with initial fluvial topography.(b) The landscape after 100 m average glacial erosion (using m = 2 in Eq. 9).Several characteristic glacial landforms are evident in the eroded landscape.Highlighted with numbered arrows are (1) flattened valley floor, (2) a hanging valley, and (3) a truncated spur.The main valley is widened forming a U-shaped cross section with steep slopes along the sides.(c) The total bedrock erosion.Erosion is concentrated in the main trunk valley along the truncated spurs and at high elevation near the headwalls of tributary glaciers.

Figure 9 .
Figure 9. Evolution of mean sliding rate as a function of mean erosion.(a) Shows results for three different sliding laws using a nonlinear erosion law (m = 2 in Eq. 9) and 80 % floatation.(b) The same as (a) but for m = 1.(c) The evolution of mean sliding for the Coulomb-friction sliding law and three different levels of floatation.The means of sliding rate erosion are computed only for the trunk valley, defined as all glaciated cells below 750 m elevation in the initial fluvial landscape.

Table 1 .
Model parameters used for all experiments.