Effect of higher-order stress gradients on the centennial mass evolution of the Greenland ice sheet

We use a three-dimensional thermo-mechanically coupled model of the Greenland ice sheet to assess the effects of marginal perturbations on volume changes on centennial timescales. The model is designed to allow for five ice dynamic formulations using different approximations to the force balance. The standard model is based on the shallow ice approximation for both ice deformation and basal sliding. A second model version relies on a higher-order Blatter/Pattyn type of core that resolves effects from gradients in longitudinal stresses and transverse horizontal shearing, i.e. membrane-like stresses. Together with three intermediate model versions, these five versions allow for gradually more dynamic feedbacks from membrane stresses. Idealised experiments are conducted on various resolutions to compare the time-dependent response to imposed accelerations at the marine ice front. If such marginal accelerations are to have an appreciable effect on total mass loss on a century timescale, a fast mechanism to transmit such perturbations inland is required. While the forcing is independent of the model version, inclusion of direct horizontal coupling allows the initial speed-up to reach several tens of kilometres inland. Within one century, effects from gradients in membrane stress alter the inland signal propagation and transmit additional dynamic thinning to the ice sheet interior. But the centennial overall volume loss differs only by some percents from the standard model, as the dominant response is a diffusive inland propagation of geometric changes. For the experiments considered, this volume response is even attenuated by direct horizontal coupling. The reason is a faster adjustment of the sliding regime by instant stress transmission in models that account for the effect of membrane stresses. Ultimately, horizontal coupling decreases the reaction time to perturbations at the ice sheet margin. These findings suggest that for modelling the mass evolution of a large-scale ice sheet, effects from diffusive geometric adjustments dominate effects from successively more complete dynamic approaches.


Introduction
The Greenland ice sheet (GrIS) is a highly dynamic system where net accumulation is compensated by both surface ablation and ice discharge at the marine margin.The system is closed by ice motion and it was essentially in balance in the decade preceding this millennium (Zwally et al., 2005;van den Broeke et al., 2009;Rignot et al., 2011).For the present state, this mass turn-over is stated to be out of balance and recent estimates for current mass loss rates based on gravimetry (Luthke et al., 2006;Velicogna et al., 2009;Schrama and Wouters, 2011), on radar/laser altimetry (Slobbe et al., 2008;Sørensen et al., 2011;Zwally et al., 2011;Moon et al., 2012), and on approaches combining discharge estimates from satellite interferometry with surface mass balance modelling (van den Broeke et al., 2009;Rignot et al., 2011) range from 0.3-0.8mm s.l.e.yr −1 (sea level equivalent) (110-290 Gt yr −1 ).About 60 % of this accelerated mass loss is attributed to increased runoff, while observed outlet glacier speed-up export the other 40 % (Rignot et al., 2011;Moon et al., 2012;Sasgen et al., 2012).Although outlet glacier thinning and retreat is observed all around Greenland (Thomas et al., 2009), the temporal and spatial pattern for the associated speed-up is far from homogeneous.For some coastal areas, observations point to regionally linked glacier accelerations (Howat et al., 2008; Published by Copernicus Publications on behalf of the European Geosciences Union.Howat and Eddy, 2011) while for other areas, a more erratic behaviour is found (Howat et al., 2010;Joughin et al., 2010;McFadden et al., 2011;Moon et al., 2012).This inhomogeneity impedes the identification of the most relevant trigger mechanisms.Hypotheses range from enhanced basal lubrication by increased summer melt (Zwally et al., 2002;Joughin et al., 2008;Bartholomew et al., 2010;Sundal et al., 2011;Schoof, 2010), to cryo-hydrologic warming by melt drainage (Phillips et al., 2010), ocean-induced grounding line retreat (Schoof, 2007;Briner et al., 2009) or enhanced calving activity (Benn et al., 2007;Nick et al., 2010).Loss of buttressing at the glacier front is believed to be induced either by local surface melt (Scambos et al., 2000), by hydrofracturing and subsequent front destabilisation, or by intrusion of warm, saline ocean water into the local fjord systems, causing submarine melting (Holland et al., 2008;Amundson et al., 2010;Straneo et al., 2010Straneo et al., , 2011;;Motyka et al., 2011;Christofferson et al., 2011).Except for cyclic surge behaviour, observations indicate that accelerations mostly originate in the marine-terminated ablation area.This area constitutes a small fraction of the entire GrIS and instant effects from speed-up remain confined to the periphery.Such marginal perturbations are directly and indirectly transmitted inland by ice dynamics and the efficiency of this signal transmission controls the impact on the interior (Nick et al., 2009).One indirect mechanism for signal transmission is diffusive surface elevation adjustment.Marginal perturbations alter the local ice thickness and surface slope, which in turn affects the flow in the immediate vicinity.The surrounding geometry consequently adjusts and thereby the signal is gradually transmitted inland.The effects of such diffusive geometry adjustments are well understood and have been studied in models for large-scale ice sheets both on long and short time-scales (Huybrechts and de Wolde, 1999;Ritz et al., 2001).However, there are also direct mechanisms that can instantaneously transmit marginal perturbations inland.One potential candidate for direct signal transmission is inherent in the force balance within a solid body of ice.Gradients in membrane-like stresses (cf.Hindmarsh, 2006) allow for such direct horizontal coupling.These stresses are often referred to as "longitudinal stresses" owing to an extensive study on flow bands.
From a theoretical point, longitudinal coupling in an ice sheet is efficient within a domain of four to ten times the ice thickness (Kamb and Echelmeyer, 1986).For the Greenland ice sheet we therefore expect values of 5 to 40 km.For a flow band setup on Jakobshavn Isbrae, Price et al. (2008) succeeded to explain instant accelerations observed at Swiss Camp by longitudinal transmission of near margin perturbations.With a linear flow and sliding law, they solve the full force balance equation and estimate that horizontal coupling has a 10-15 km area of direct influence.For a similar 2-D application on Helheim, Nick et al. (2009) highlight the transient propagation of perturbations initiated at the glacier terminus.The immediate response confirms a 10 km coupling radius to the applied marginal forcing.The effective longitudinal coupling length increases when accounting for the nonlinear character of ice creep and basal sliding (Kamb and Echelmeyer, 1986;Price et al., 2008).For a typical fast flowing Antarctic ice stream, it can reach up to 40 km (Williams et al., 2012).Flow band geometries are however limited to the effect of longitudinal stress gradients and omit effects from transverse shearing, which occur widely in many stream-like outlet glaciers embedded in the GrIS.Moreover, the purpose of such flow line setups is often to reproduce observed flow evolution rather than the study of the influence of direct horizontal coupling on the geometry evolution.
In this paper, we apply a 3-D thermo-mechanically coupled Greenland ice flow model that includes the option to activate direct stress transmission via gradients in longitudinal stresses and in transverse horizontal shearing.First, the model is described in Sect. 2 together with a justification for its applicability to direct stress transmission.Section 3 sets the experimental frame for the centennial perturbation scenarios on Greenland, while results are subsequently discussed in Sect. 4. Introduction of an instant reaction timescale provides a tool to summarise results from all experiments for different resolutions and model versions in Sect. 5. Before the final conclusion, a discussion of the entire model setup assesses the reliability of the results.

Stress terminology
The fundamental force balance is presented in Appendix A and serves as a basis for the used stress terminology.Based on scaling arguments, vertical plane shearing (σ xz , σ yz ) exerts major control on the dynamics of large-scale grounded ice bodies, as exemplified in the shallow ice approximation (SIA; Hutter, 1983).Following the terminology in Hindmarsh (2006), normal stresses acting in the horizontal direction (σ xx , σ yy ) and horizontal shearing on vertical planes (σ xy , σ yx ) are referred to as membrane-like or simply membrane stresses.Here we refer to σ xx , σ yy as the longitudinal stresses acting in the horizontal direction parallel to the normal vector of the respectively addressed surface.The other two membrane stresses we refer to as transverse horizontal shear stresses (σ xy , σ yx ), or simply transverse horizontal shearing.Their influence on ice dynamics is referred to as direct horizontal coupling or direct stress/signal transmission.Note that this terminology is linked to the imposed coordinate system and is best applicable to large-scale geometries.For generalisation to more complex geometries, a consistent and precise terminology will depend on local geometry.

Abbreviation Experiment description
MarAsl2 Doubling the sliding factor upstream of marine terminated outlets.TotAsl2 Doubling the sliding factor over the entire sliding area.MarCut Singular calving event for a marine margin extending 20 km inland.

General model characteristics
The three-dimensional thermo-mechanically coupled ice sheet model comprises three main components that respectively describe the ice rheology, the Earth crust's isostatic rebound and the mass balance at the upper and lower ice sheet boundaries (Huybrechts and de Wolde, 1999;Huybrechts, 2002; to which the interested reader is referred for a detailed description of all model components).The optional higher-order (HO) ice-dynamic core is described in Fürst et al. (2011).The surface mass balance model is based on the widely used positive degree-day/retention method.As input serves a generic temperature field, which depends on latitude and surface elevation, together with a temperature scaled precipitation field (Huybrechts, 2002).Ice dynamics, isostatic adjustment and mass balance are implemented on a computational grid with 5 km horizontal resolution and 30 nonequidistant layers.The vertical grid spacing is refined towards the bottom where vertical plane shearing is concentrated.Geometric input has been modified from the original Bamber et al. (2001) data as described in Huybrechts et al. (2011).

Model versions for ice dynamics
The ice-dynamic model component is based on a viscous rheology assuming isotropy within the ice body.Ice flow results from internal deformation together with sliding over soft till or bedrock, but is limited to areas where basal temperatures are within 1 • C of the pressure melting point.Modelled ice creep follows a Glen-type constitutive relation (Eq.A3) with exponent three linking the stress field to strain rates and velocities (parameters in Table 1).The used proportionality factor depends on ice temperature and ice age.Together with the force balance (Eq.A1) and the incompressibility Eq. (A2), they form the most general set of equations for ice deformation, solved in so-called Full-Stokes (FS) models (Zwinger et al., 2007;Jouvet et al., 2009;Seddik et al., 2012).In our approach, this system of equations is not solved in its full complexity but different approximations are applied concerning ice deformation and basal resistance.This results in altogether 5 different model versions (cf.For ice deformation our model optionally applies the classical SIA approach or an LMLa higher-order (HO) core (Fürst et al., 2011).The latter version captures dynamic feedbacks arising from gradients in membrane stresses (cf.Sect.2.1).These two options are combined with either of three possible descriptions of basal resistance.A first approach assumes that basal resistance locally supports the driving stress (DR).A second approach uses the shallow shelf approximation as a sliding law (ME) (cf.Ritz et al., 2001;Bueler and Brown, 2009;Schoof and Hindmarsh, 2010).In this case, resistance at the base is dominated by membrane stresses while any effect from vertical shearing on sliding is neglected.Reconciling these two extreme approaches, a simplified resistance (SR) equation includes effects from both vertical plane shearing and membrane stresses at the base.The simplification used in this third approach is inherited from the HO ice deformation and stems from a scale analysis of gradients in the velocity field (cf.Appendix B).
Table 2 summarizes the five model versions that arise from combining these approximations used for ice creep and basal resistance.Note that the SR SIA model version actually needs a full determination of the 3-D higher-order solution in order to extract the velocity field at the base.In addition, the ME HO model version is not considered here because membrane stress activity for the basal and deformational part of ice dynamics is more self-consistent in SR HO.

Experimental design
The experimental frame for this study is set by three schematic perturbation scenarios and focuses on the centennial response of the GrIS.These experiments serve mainly as an intercomparison study for the five model versions with different dynamic complexity but they also give indications on the dynamic ice loss within one century.Optional grid spacing of 20, 10 or 5 km allows us to investigate the impact of resolution on capturing direct signal transmission.

Initialisation
An interglacial equilibrium state serves as the standard initial condition.The climatic forcing is obtained from a generic temperature field that depends on latitude and surface elevation (details in Huybrechts, 2002).The initialisation on 5 km takes 25 kyr and is preceded by a 200 and a 100 kyr equilibrium run on respectively 20 and 10 km.For this spin-up the DR SIA model setup is chosen and we denote this standard interglacial with "IS".For assessing the robustness of the results under different initial conditions, we experimented with two other spin-up methods.The first one uses the SR HO model version to produce the interglacial equilibrium.It starts from the respective IS geometry and is evolved for either 50, 20 and 5 kyr to an initial state "IH".For the third initialisation method, the model is spun up over the last two glacial cycles in the same way as described in Huybrechts et al. (2011).The resolution was successively increased from 20 km to 10 km at 20 kyr BP and to 5 km at 3 kyr BP.

Perturbation experiments
We conduct three experiments that address various aspects of Greenland's ice dynamic response (cf.Table 3).Their focus is on the century timescale while they are designed in a way that each model version is applicable.The central experiment prescribes a step function in the basal sliding coefficient A sl but spatially confines this amplification to the marine terminated boundary (experiment MarAsl2).Our choice for an A sl amplification factor of 2 in grid boxes within 40 km upstream of the marine boundary is loosely based on the details of the 1992/2002 speed-up of Jakobshavn Isbrae (Joughin et al., 2008).The dynamic complexity of our model versions is expected to influence the inland propagation of geometric adjustments to such a marginal perturbation.In a second experiment (TotAsl2), the basal sliding coefficient is again doubled but now over the entire sliding area.This experiment can be considered as an extreme lubrication event, which by far exceeds the ablation area.In a last experiment, the model response to a massive singular ice removal event at the marine margin around Greenland is investigated.Grid cells that have a marine boundary are instantaneously emptied.In order to ensure comparability between all resolutions ice removal is enforced up to 20 km inland.
To enable a clean comparison between the five model versions, ice velocities in the first grounded point are based on the DR SIA model version.This provides unified perturbations for all experiments.Otherwise we would introduce numerical differences in the implementation of the lateral boundary conditions originating from the non-local influence on the velocity solution in model versions with horizontal coupling, which requires zero velocity outside the ice sheet.Although this approach inhibits membrane stress activity in the first grounded ice point, none of the used resolutions actually resolve the complex stress situation in these zones.Additional tests (not shown) brought to light that applying the DR SIA lateral boundary condition at the first grounded point either all around Greenland or restricted to the marine margin has no significant influence on the presented results.Using SR HO at the lateral boundary, the qualitative differences in volume evolution and in thickness response remain the same between models that use DR SIA and SR HO up-stream.However, the centennial mass loss increases by 17 % with this second boundary condition but independently of model version.Unperturbed control runs are conducted to assess the effects of model adjustment and background evolution.

Results
Since the results were qualitatively similar for each initialisation technique, only results from the IS spin-up are presented in the following.The interglacial ice sheet extent and volume are similar to the present-day configuration (Fig. 1).A notable mismatch is the excessive amount of ice on Peary Land because bed troughs are not well resolved and because of shortcomings in the surface mass balance field during the spin-up procedure.The associated velocity field on 5 km resolution reproduces many narrow stream features that appear at locations where recent observations show equivalent fast ice flow (cf.Fig. 1 and Joughin et al., 2010).For the 20 km initial state most of these narrow streams are absent, which www.the-cryosphere.net/7/183/2013/The Cryosphere, 7, 183-199, 2013 underlines the gain from grid refinement.Note that the modelled stream features arise naturally from resolved details in bedrock topography and subglacial heating.They are not the result from optimising the basal sliding coefficient to match the observed velocity field as in Price et al. (2011).
In our setup, the IS surface velocity field is in equilibrium with the DR SIA equations and the continuity equation.The velocity field computed with the SR HO variant for the same geometry largely resembles the DR SIA solution.Effects from higher-order dynamics only become apparent in minor details of the ice flow (cf.Fig. 1).The DR SIA velocity magnitude exhibits sharper jagged features that are effectively smoothed by direct horizontal coupling in the full higher-order version.Moreover the ice divide velocities tend to be higher in the SR HO version as a result of longitudinal stress transmission.Furthermore, effects from transverse horizontal shear are readily visible across narrow ice streams, where fast stream velocities couple with the stagnant ice in the surroundings (see e.g.Kangerlussauq, Fig. 1c, d).
In the following the transient effect of direct horizontal coupling is studied on various resolutions.First the GrIS volume evolution is determined for all perturbation experiments (cf.Table 3).The GrIS mass loss for MarAsl2 is then decomposed in a purely dynamic signal and a component arising from the feedback with surface mass balance.The dominant dynamic feedback is further investigated by comparing all model versions.Finally, results covering all model versions, grid resolutions and marginal perturbation experiments are analysed by introducing a reaction timescale to characterise their individual response behaviour.

Dynamic mass loss
For all experiments, the total volume loss of the GrIS after one century falls into a range of 1 and 12 cm eustatic sealevel equivalent (s.l.e.), cf.Fig. 2.This centimetre volume response is in line with estimates for present-day rates of increased ice export of about 0.2-0.4mm s.l.e.yr −1 (van den Broeke et al., 2009) and projections for future dynamic response (Price et al., 2011).The results have to be interpreted relative to Greenland's total mass balance.The modelled net accumulation for the IS equilibrium geometry is equivalent to 1.5 mm s.l.e.yr −1 .It is balanced by surface runoff and ice discharge at rates of 0.84 (55 %) and 0.66 mm s.l.e.yr −1 (45 %), respectively.Doubling of the sliding coefficient in DR SIA initially almost doubles the ice discharge and would therefore extrapolate to a centennial mass loss of about 66 mm s.l.e. in case such rates could be sustained.However, all experiments show that initial mass loss rates are not sustained.Even for a doubling of sliding over the entire sliding area (TotAsl2), upstream supply of ice is not sufficient to sustain the discharge rates.They gradually decrease and add up to a mass loss of 40 mm s.l.e. after one century (Fig. 2c).For MarAsl2 (Fig. 2a), initial mass loss rates abate even faster due to a lack of upstream inflow, adding up to 16 mm after one century.Similarly, the initial release of some 40 mm s.l.e.ice volume in the MarCut experiment is followed by a quick decrease in mass loss rates (Fig. 2d).Rerunning the MarAsl2 experiment from the IH interglacial or from the present PS state (cf.Fig. 2b) shows a qualitatively similar mass evolution.On 20 km, the MarAsl2 experiment was repeated with amplification factors of five and ten (results not shown).The results indicate that the centennial mass loss scales almost linearly for such forcing magnitudes.
A comparison of the different model versions shows that the triggered mass losses are similar (Fig. 2).Differences lie within 20 % of the total signal, so that additional effects from the ice-dynamic complexity of each model version are generally small.For all experiments, an increase in dynamic complexity actually reduces the centennial mass loss.This negative feedback is especially strong when activating a full nonlinear sliding law (dashed lines in Fig. 2a) reaching 20 % of the total signal.This highlights the importance of sliding in attenuating the mass loss signal.This mass loss reduction is reproduced on all resolutions (cf.Sect.4.4) and is robust under the various initial states.
The fact that the centennial mass loss reduces as the ice dynamic complexity increases raises the question whether this is caused by ice dynamics alone or if there is a dominant feedback with surface mass balance.Differences in ice dynamics could increase thinning rates and together with the elevation feedback increase surface runoff (Huybrechts et al., 2002).In the MarAsl2 experiment, this feedback explains at most 10 % of the mass loss differences between our model versions (cf.Figs. 3 and 2a).Differences in this feedback between model versions are not consistent neither over the various resolutions nor when allowing for direct horizontal coupling.Therefore it can be excluded from being the decisive factor here.Differences in mass loss are rather the consequence of reduced discharge at the marine ice front.

Inland propagation of marginal perturbations
For MarAsl2, the modelled reduction of ice discharge at the marine margin depends on the details of the upstream propagation of these accelerations.To investigate spatial differences in successive geometric adjustment, we first analyse the initial velocity response.We conclude this section with a decomposition of the velocity response into its basal and deformational component.
The basal sliding coefficient A sl is doubled within an area extending 40 km inland from the marine ice front.The initial response of the velocity field under identical geometries is shown in Fig. 4 for the three main and more commonly used model versions.In the DR SIA setup (Fig. 4a, d, e), accelerations are obviously limited to the actually perturbed grid points and allow a delineation of the A sl amplification zone.Note that the perturbation is additionally confined to the actual extent of the sliding area, which becomes visible for the close-ups on Jakobshavn Isbrae (Fig. 4d) and Kangerlussauq (Fig. 4e).Including horizontal coupling, ME SIA (Fig. 4b, f, g) and SR HO (Fig. 4c, h, i) reveal how the initial acceleration instantaneously radiates some distance inland.Transverse horizontal shearing affects a vicinity of about 5-10 km for both model versions, whereas the cou-pling length for longitudinal stress gradients seems larger.ME SIA exhibits values of about 20 km while SR HO expands the coupling to more than 40 km (cf.Fig. 4f, g, h, i).For ice thicknesses of about one kilometre, these values are in broad agreement with the theoretical coupling length (Kamb and Echelmeyer, 1986).A greater coupling length is possible for high sliding or elongated, deep bed troughs as is the case for Kangerlussauq (Fig. 4e, g, i).For Jakobshavn Isbrae, Price et al. (2008) exemplarily demonstrate the positive influence of basal sliding on the longitudinal coupling length.
Figure 5 illustrates spatial differences in inland propagation of thickness changes.For the DR SIA model, geometric changes are a result of simple diffusive thinning (upper panels in Fig. 5).Thickness changes in MarAsl2 are initiated at the marine margin and are transmitted up to the divide, where elevation only changes by some millimetres after one century.In the first ten years, diffusive adjustment can locally result in peripheral thickening by intense upstream inflow (upper panels Fig. 5).However in a later stage, thinning along the marine margins is predominant while sporadic thickening is observed locally, mainly in the east.Such thickening stems from complicated bed and flow geometries of individual outlet glaciers.Since diffusive thinning dominates the transient response in all model versions, their thickness evolutions are given with respect to the DR SIA reference.Initially, direct horizontal coupling hardly affects the geometric adjustment beyond the area of diffusive thinning.But close to the marine margin horizontal coupling causes more pronounced thinning.In a later stage, its integrated effect shows up as an extra thinning signal in the interior, most pronounced for model versions that resolve vertical variations in membrane stresses (DR HO and SR HO).Yet in all models, a general thickening pattern emerges in the upper regions of the ablation zone, in between the zones of relative thinning.This spatial pattern appears to be a fingerprint of accounting for horizontal gradients in membrane stresses, especially when resolving them vertically (SR HO and DR HO).Though not decisively altering the mass loss response on centennial timescales, the additional thinning reduces the ice export and thus provides an explanation for the decreased total mass loss (cf.Fig. 2).The relative thickening of the upper ablation areas arises from downstream damming by the thinned margin.But also upstream accelerations (see Fig. 6) add to this ice bulge, in particular for SR HO and DR HO.Since this relative thickening is situated in the upper ablation region, a positive feedback between surface elevation and mass balance amplifies this phenomenon (cf.also Huybrechts et al., 2002).After 50 yr of geometry evolution, our model versions branch out significantly and velocity fields also hold information on differences in integrated ice thickness changes (cf.Fig. 6).Prominent are the relative deceleration in the region of dynamic thickening and the marginal speed-up.ME SIA restricts the effect from membrane stress gradients to the basal layer and velocity differences are dictated by the basal component.Consequently the deformational component hardly exceeds the area of basal velocity changes (Fig. 6a, c) and additional thinning further upstream remains negligible.Interior thinning however is experienced in the SR HO model, where direct horizontal coupling throughout the entire vertical ice column transmits the acceleration beyond the sliding area extent (see Fig. 6b, d).
Altogether, the spatial propagation of marginal perturbations confirms that dynamic complexity indeed alters the geometry evolution, but this is far outweighed by the dominant background thickness diffusion.Significant effects on geometric adjustment from direct horizontal stress transmission cannot be confirmed even in the case of full nonlinear sliding.

Decomposition of mass loss attenuation
Direct horizontal coupling causes additional thinning along Greenland's marine margins for the MarAsl2 perturbation.Whether this thinning alone is sufficient to explain the observed mass loss differences between model versions remains yet unclear.For this purpose, we decompose the ice discharge into a contribution from thinning alone and from acceleration alone.Since discharge is proportional to the ice thickness H and the velocity magnitude V near the marine terminated margin, the ice discharge difference D between a forced and a reference experiment becomes: The discharge difference D is a cumulative quantity that sums up differences between products of a mean ice thickness H and a mean velocity V .The spatial mean is found by covering all the forced ice points while velocities are vertically averaged.The proportionality factor C has the units of a length scale and is proportional to the chosen grid spacing.The discharge difference between a perturbed and a reference run consequently depends on three components: differences in ice thickness evolution H = HP (t) − HR (t), differences in velocity evolution V , and a combination of both.Their respective signs, shown in Fig. 7, confirm an average increase in velocity magnitude ( V > 0) and thinning along the margin ( H < 0).In addition, their magnitudes indicate the dominance of velocity differences for the ice discharge evolution.The cumulative effect from differences in ice thickness evolution alone is one order of magnitude lower than the effect from velocity changes.Accounting for membrane stresses in ice dynamics amplifies the marginal thin-ning (cf.Fig. 7a) and therefore partially explains the negative feedback on the centennial mass loss (cf.Fig. 2).The combined effect from changes in both velocity and ice thickness (cf.Fig. 7c) adds a small and negative contribution to ice discharge.But the mass loss in different model versions ultimately diverges because of a different velocity evolution (Fig. 7b).Since the marine margin has local DR SIA velocities, this effect must arise from upstream, where different flow adjustments and resultant geometric changes (with respect to DR SIA) reduce the inflow and thus cause extra thinning and deceleration.
For a further analysis of the MarAsl2 velocity response with emphasis on upstream ice supply, we focus on the sliding area of Greenland below the 2000 m topographic contour.The resulting region provides the main supply for ice discharge at the marine margin.We additionally decompose the velocity field into its sliding component and the vertically averaged, deformational component (Fig. 8).Both components show an instant increase after the perturbation is applied, which on average is about 250 m yr −1 for sliding and 17 m yr −1 for deformation.Yet the acceleration is attenuated fast in both components.The deformational component almost completely returns to the reference value within one hundred years, while sliding remains at about 120 m yr −1 above the reference.Already the difference in magnitude highlights the importance of the sliding component for horizontal coupling.For the attenuation of the initial basal acceleration (Fig. 8a), it is striking that all model versions with horizontal coupling at the base (ME SIA, SR SIA, SR HO) group and show a more expressed dampening of the initial perturbation.An integration of this sliding perturbation over time would explain most of the model differences seen in the total mass loss (Fig. 2a).It seems only logical that for a perturbation in the sliding coefficient, it is the sliding velocity that controls the response.Yet for the geometry perturbation in the MarCut experiment, this dominance is confirmed.It presumably results from our choice of the respective parameters in the sliding law (Eq.B3) and the flow (Eq.A3).Also DR SIA and DR HO follow a similar sliding attenuation, which is not the case for their deformational component (Fig. 8b).Here model versions group according to the used deformational approximation with a weaker dampening for HO models.In essence, differences in ice deformation near the margin neither explain the spread in cumulative discharge (Fig. 2a) nor its magnitude.Dampened ice export thus results from faster attenuation of basal accelerations near the margin when including dynamic effects from horizontal gradients in membrane stresses.

Reaction times
Accounting for the dynamic effects of membrane stresses and the resulting non-local horizontal coupling allows for a faster attenuation of the initial perturbation.To quantify this, we aim to distil differences in transient responses for each model version, each grid resolution, and each experiment into one scalar.We introduce a reaction time based on an exponential decay: The cumulative ice discharge D has a saturation value S with a mean characteristic reaction time τ R .The retrieved values are averaged over one century.In this way information from the entire mass evolution enters the two parameters.The resultant mean reaction time shows magnitudes of several decades for marginal forcing (cf.Fig. 9) but it exceeds hundreds of years for the TotAsl2 perturbation experiment.Therefore, we focus on the marginal forcing experiments with clear attenuation within our centennial scope.
A common characteristic for all margin perturbation experiments is the reduction of the reaction time when direct horizontal coupling is included in the basal layer (Fig. 9).The reason is that the instant velocity response exceeds the area of direct perturbation and thus allows geometric adjustment in a wider area than in DR SIA, resulting in faster attenuation.However, differences due to the model versions remain within 10 % of the DR SIA reference.Reaction time differences between the model with a linearised effect of membrane stress gradients in the sliding law and the one accounting for their full nonlinear impact (empty circles, Fig. 9a) are of the same order as those between other model versions.Therefore their linearisation (cf.Appendix B) is acceptable when the interest is in the centennial volume response.Yet the response behaviour is sensitive to the chosen grid resolution and reaction times vary by up to 40 %.For both marginal perturbation experiments, reducing the grid resolution lowers the reaction time.This would be an undesirable effect if it simply arose from pure grid size reduction.For this reason we repeated all experiments by simply supersampling the 20 km bed topography for the spin-up without using the detailed information from a higher resolution.Then the MarAsl2 set of experiments (not shown) has reaction times of about 40 yr with no monotone dependence on grid resolution.Consequently, the reduced reaction time in our experiments arises from additional information in the input data, which interfere with the spin-up and in turn affect the experimental setup.This behaviour is indeed a wanted effect because in reality geometries are complex and will affect the inland propagation of marginal perturbations.The MarCut experiment however shows the grid dependent reduction of the reaction time in both setups.Here this is a relict from an experimental setup disturbed by grid spacing.While MarAsl2 prescribes a constant 40 km extent of the sliding factor increase, the Mar-Cut perturbation has a direct effect only on the first grounded grid point upstream of the calving points.In summary, a consistent picture emerges for all marginal forcing experiments and for each resolution.Ice sheet reaction times are reduced when direct horizontal coupling operates and therefore allow for a faster attenuation of marginal perturbations.Yet on a century timescale, this faster attenuation does not play a large role for the total mass loss of the GrIS.

Discussion
To assess the effect of direct signal transmission or horizontal coupling on the dynamic response of the GrIS over the next century, five model versions of a 3-D ice-dynamic model were used with different dynamic complexity that span the gap between the SIA and a higher-order LMLa model for ice evolution.Although effects from membrane stress gradients and their vertical variations are captured, its application is limited by two simplifications (Pattyn, 2003).The decisive one is that the glaciostatic approximation neglects bridging effects in the vertical force balance.In addition, horizontal gradients in the vertical velocity field are assumed to be negligible compared to vertical gradients in horizontal velocities.One might question whether these two assumptions result in a substantial loss in dynamic complexity that might inhibit important feedbacks for the transient ice sheet behaviour.Based on a scaling from an asymptotic analysis, Schoof and Hindmarsh (2010) present arguments that affirm the applicability of an LMLa model for large aspect ratios.The horizontal component of the LMLa velocity solution shows an error of second order in this ratio.Since the GrIS extent is at least two orders of magnitude larger than its thickness and since the horizontal grid spacing exceeds the ice thickness, our large-scale ice sheet model remains an appropriate tool for ice dynamics.Although Schoof and Hindmarsh (2010) emphasise the validity of the LMLa approach for slip and non-slip areas separately, the general scaling arguments might be violated in the transition zone (Schoof, 2006).However, the higher-order LMLa approach is valid for dynamics in the vicinity of the ice divide where the nonlinearity of the viscous material causes a unique internal layer architecture, the Raymond effect (Raymond, 1983).In general, an LMLa model seems a feasible tool for transient, large-scale applications and it is assumed that additional feedbacks in the FS equations are limited to smaller scales (< 1 km) when stress regimes are well resolved.
Another issue concerns the poorly constrained sliding mechanisms beneath large-scale ice sheets (Fowler, 2010).Since the focus is on inland signal propagation, the applied sliding law should be capable of realistically transmitting margin perturbations upstream.Observations on Rutford Ice Stream (Gudmundsson, 2011) suggest that a nonlinear mechanism is a prerequisite to explain the upstream ice response to tidal forcing of different frequencies.A Weertman sliding law with an exponent of three reasonably explains the observed response pattern.This justifies its use for assessing the influence of direct horizontal coupling on inland signal transmission in the presented experimental framework.
At a numerical level, our choice for three distinct grid resolutions alters and potentially inhibits direct horizontal coupling.Since we do not find a strong grid dependence of the total mass loss, it remains possible that direct horizontal coupling is not resolved properly on any used resolution.However, on a 5 km grid, the initial speed-up in a 40 km vicinity (Fig. 4) is resolved and in agreement with theoretical estimates (Kamb and Echelmeyer, 1986).We are therefore convinced that membrane stress activity is captured in an acceptable way in our large-scale ice flow model.Except for resolving successively more geometric details, we do not expect that further grid refinement would drastically alter model differences in the spatially integrated mass loss for these experiments.Underlining the consistency of our results, all of them are robust under different initial conditions (Fig. 2b).
A last point concerns the artificial nature of our perturbation experiments.All of them cause a dynamically dominated response within one century.But the initial perturbations have in each case a rather extreme character.The MarAsl2 experiment mimics a simultaneous speed-up of the marginal zones of all outlet glaciers around the GrIS while observations indicate a more variable behaviour (Moon et al., 2012).In a similar way, the TotAsl2 experiment doubles velocities within the entire sliding area, as an attempt to simulate enhanced lubrication due to increased surface melt.The  2).This timescale reflects the immediate response to a stepwise perturbation rather than a full ice sheet equilibration of many thousands of years.spatial extent, the duration and the magnitude of the amplification in TotAsl2 exceed recent observations by far (McFadden et al., 2011;Sundal et al., 2011).A simultaneous speedup is also imposed for MarCut.From the three experiments this is the most severe one, with an equivalent sea level contribution that exceeds 100 mm.The other two experiments result in 45 mm and 16 mm after 100 yr.Though each experimental perturbation is physically extreme, the resulting sea level contribution is in a realistic range with respect to present estimates for the increase in ice discharge rates (van den Broeke et al., 2009;Rignot et al., 2011).The reason is a fast attenuation of the initial perturbation and the consequently reduced centennial response.Therefore these experiments should be judged from the projected volume response rather than from the applied physical perturbation.The modelled mass loss range also agrees with other studies on dynamic response of the GrIS (Price et al., 2011;Graversen et al., 2010).Graversen et al. (2010) conducted future projections under various climate scenarios, parameterising a doubling of sliding speeds around Greenland.The resulting mass loss in these projections does not exceed 17 cm by the end of the 21st century, including uncertainties from the climate model ensemble.In their study, 15 % or 25 mm of the total mass loss are associated with outlet glacier accelerations, which is somewhat lower than in TotAsl2.The reason for this is the imposed future warming that causes additional surface melt near the margin, where it removes the ice before it can reach the marine margin for calving.

Conclusion
A thermo-mechanically coupled 3-D large-scale ice sheet model with five versions of dynamic complexity was used to assess the influence of membrane stress gradients on the centennial response of the GrIS under various marginal perturbations.Results indicate that the SIA, which does not include effects from direct horizontal coupling, already captures the major response characteristics.In fact, the inclusion of membrane stress gradients actually reduces the volume loss by up to 20 %.The reason is a faster attenuation of these perturbations by instant upstream accelerations and successive geometric adjustments.Including direct horizontal coupling also increases the area of direct influence.Our results exclude a dominant feedback between surface elevation changes and surface mass balance; it is rather the ice discharge that abates more quickly and prevents sustained mass loss rates.The discharge reduction is moreover not explicable by additional marginal thinning but rather by faster attenuation of the speed-up signal.Furthermore, the experiments have identified sliding as the component that reacts fastest and explains the differences in mass evolution.In summary, non-local effects from membrane stresses allow for a faster readjustment and attenuation of margin perturbations.Yet the general response of all model versions is similar to DR SIA, where signal transmission is accomplished by thickness diffusion alone.
An analysis of all experiments, model versions and grid resolutions shows that the reaction timescale is primarily dependent on the applied perturbation.The stronger the perturbation and the larger its spatial extent, the longer the reaction time.Putting aside the resolution dependence, the reaction time decreases by not more than 10 % when including direct horizontal coupling.This underlines the small effect of the dynamic complexity on the centennial volume evolution of the GrIS.
Our findings have important consequences to decide on the degree of complexity needed when studying the dynamic response to prescribed marginal perturbations.If the interest is in the volume evolution of an ice sheet, the classical SIA approach is sufficient to capture the main characteristics of the inland response.The higher-order approach SR HO is dynamically more complete but the additional computational cost (factor very dependent on chosen accuracies but in our case about 100) is not warranted by the limited gain in ice sheet volume response.A compromise could be the ME SIA approach (Bueler and Brown, 2009) with direct horizontal coupling at the base only.This approach is computationally much more efficient than SR HO (factor 20).Since basal velocity adjustment controls the signal attenuation ME SIA seems more appropriate than DR HO, though computationally they are comparable.When SR HO or FS models become computationally too expensive for a specific setup, we recommend the use of an ME SIA equivalent approach (Hindmarsh, 2006;Bueler and Brown, 2009;Schoof and Hindmarsh, 2010).
In summary, differences in inland propagation of margin perturbations due to the degree of simplification of the force balance affect the ice sheet volume response only marginally.A self-accelerating dynamic feedback that would facilitate interior drainage and ultimately could trigger ice sheet disintegration is not confirmed.In addition, all perturbation experiments have a rather extreme character but the resulting mass loss is very limited.Simultaneous speed doubling all along the marine margin or in the entire sliding area goes far beyond observed changes during the last decade.This strengthens the view that the dominant effect on volume change of the Greenland ice sheet for the 21st century will come from surface mass balance changes.The largest source of uncertainty for future projections of ice sheets therefore arises from the climatic projections used as input rather than from the complexity of the ice-dynamic model.

Variants for ice deformation and basal resistance A1 Ice deformation
Our glacial system is placed into an orthogonal coordinate system with three unit vectors {e x , e y , e z } in respectively horizontal x-, y-and vertical z-direction.The vertical axis of our coordinate system is chosen to be perpendicular to isolines of the gravitational field.In essence ice deformation is based on two balance equations for mass and momentum combined with a constitutive relation that links the stress tensor σ to strain rates and ultimately to the velocity field u = (u x , u y , u z ).Considering an incompressible ice body, the balance equations take the following form: Here, the gravitational acceleration g and the ice density ρ are assumed constant (respective values are provided in Table 1).Ice creep is controlled by deviations from cryostatic pressure (τ kl = σ kl − 1 3 δ kl m σ mm with k, l, m ∈ {x, y, z}) rather than by the mean pressure itself.A simple power law is used to relate deviatoric stresses to strain rates ε for polycrystalline ice. ) This constitutive equation is often referred to as Glen's flow law while the effective viscosity η comprises the nonlinear dependence in the system because it itself depends on the effective strain rate εe , a combination of all strain components εkl = 1 2 (∂ l v k + ∂ k v l ).The dependence on pressure melting point-corrected temperature T • of the rate factor follows an Arrhenius type function A (T • ) = A 1/2 exp(−B 1/2 /(RT • )) using two parameter sets for different temperature regimes.The rate factor also shows an age dependence since Wisconsin ice is assumed to deform more readily than Holocene ice for similar stress and temperature conditions.This rheologic influence of the ice age is described in Huybrechts (1994).
Our model offers the possibility to choose between two approximations for solving the Full Stokes (FS) equations (Eqs.A1, A2 and A3).The first is the shallow ice approximation (SIA), which is standard practice in modelling largescale ice sheets.In this approximation any influence from membrane stresses (as defined in Sect.2.1) is neglected, which results in a balance between driving stress and vertical plane shearing.
Here s indicates the surface elevation.Optionally we can determine the stress and velocity field by switching to a higher-order variant (HO) specified as LMLa (Hindmarsh, 2004;Schoof and Hindmarsh, 2010;Dukowicz et al., 2010).In this approximation of the force balance, horizontal variations in lateral vertical shearing are negligible in the vertical balance Eq.(A1).In other words, the ice overburden pressure is locally compensated, meaning that any bridging effects supporting ice weight sideways are comparably small.
www.the-cryosphere.net/7/183/2013/The Cryosphere, 7, 183-199, 2013 This variant of a higher-order model holds a second simplification that makes the force balance independent of vertical velocities and thus allows to solely solve for the horizontal components.This decoupling becomes possible by considering horizontal gradients in the vertical velocity field, small compared to vertical gradients in the horizontal components ∂ i v z ∂ z v i (cf.Eq.A8).The vertical velocity component is in turn determined by mass conservation, i.e. incompressibility (Eq.A2).As lateral boundary conditions to this force balance, ice thicknesses and velocity components are set to zero.

A2 Basal resistance
Ice deformation from SIA and HO are completed by three options to determine the basal resistance (see Table 2).In general the basal resistance τ bi compensates effects from vertical plane shearing and stress transmitted by membrane stresses.
Here sub-and superscript b indicate that stress terms have to be evaluated at bed elevation (z = b).τ || bi is the stress component parallel to the bed surface.A first option for basal resistance arises from assumptions made for the SIA.Neglecting effects from membrane stresses, basal drag is balanced by vertical plane shearing at the base, which equals the driving stress τ d .
Calculating the basal drag in this way is referred to as the driving stress (DR) option.In contrast to this dominance of vertical plane shearing at the base is an approach that highlights the importance of membrane stresses.It stems from assuming a vertically uniform stress and velocity field, thus without vertical plane shearing, while the whole setup is balanced by membrane stresses.In analogy to Bueler and Brown (2009) and using identical frictional forcing at the ice surface and at its base, the resulting force balance is used to determine the basal stress field.This resistance option with dominant membrane stresses (ME) resembles a shallow shelf approximation (SSA) apart that for our purpose it serves to determine the condition for the frictional basal boundary.As a result from vertical integration the effective viscosity η SSA is purely determined by horizontal strains that control the evolution of a floating shelf.In analogy to the HO ice deformation, a last option arises from neglecting horizontal gradients in the vertical velocity field for the strain rates.As a consequence a HO version of the resistance Eq. (A9) is obtained.This way of balancing the basal drag uses a slightly simplified resistance (SR) equation and is referred to as SR option.Together with two options for ice deformation, this gives six versions of treating basal resistance and deformation in our model (see Table 2).

Boundary conditions
In order to find a unique solution for Eqs.(A1), (A2) and (A3), boundary conditions are required.At the ice-free points around the lateral boundary we set not only the ice thickness H to zero but also the velocity field.This Dirichlet boundary condition is widely used in ice sheet modelling although resulting margin gradients are in essence dictated by grid spacing.For the vertical boundaries, conditions become more elaborate.

B1 Free surface
The upper surface boundary condition is set up in a similar way as the one for the base Eq.(A9) except that it must be evaluated at surface elevation (z = s).for i = j and i, j ∈ {x, y} However, since the contact with the atmosphere does not provide any resistance to the ice flow, its drag is zero, imposing an internal balance of stresses at the ice surface.
Note that for the SIA this relation is automatically fulfilled since vertical shearing becomes zero at the upper boundary.

B2 Basal sliding
In order to close the system of equations, boundary values have to be prescribed for the basal drag in Eq. (A9).In other words, one has to find a relation for basal sliding dependent on basal characteristics, whose influence is elaborate and still controversial (cf.Fowler, 2010).A main distinction has to be made between a pure ice bedrock contact and a case where the ice rests on a layer of till.In the first case, sliding is controlled by regelation, enhanced creep, effective water pressure, bedrock roughness and cavitation geometry (Fowler, 2010;Schoof, 2005Schoof, , 2010)).When ice overlays a till layer, finding an appropriate sliding relation becomes more elaborate since a reasonable rheology for till must be formulated.As a granular material with a yield stress, observations confirm that till is well described assuming a plastic rheology (e.g.Tulaczyk, 2000).However, field measurements on Iceland (Boulton and Hindmarsh, 1987) indicate a nonlinear viscous rheology with an exponent of three (m = 3).Such a nonlinear mechanism was found to be a prerequisite for realistic inland transmission of tidal forcing in Antarctic ice streams (Gudmundsson, 2011).
A sl H τ m−1 b τ bi with i ∈ {x, y} (B3) Here u b is the basal velocity magnitude and u bi each vector component.In the literature, this viscous approach is expressed either in a nonlinear basal friction equation or a nonlinear sliding relation.C is the friction parameter, while β 2 is the friction coefficient and both are inversely related to the basal sliding coefficient A sl .An enhancement of sliding towards the ice sheet margin is introduced by a weighting of the coefficients with ice thickness.Equation (B3) together with Eq. (A9) allow for a solution of the partial differential Eqs.(A1), (A2) and (A3) uniquely.However for practical reasons a partially linearised sliding relation is introduced.
Here the nonlinearity of the sliding relation is accounted for by the driving stress τ d , which is fully prescribed by the ice geometry.This leaves a linear relation between basal velocities u b and basal drag τ b with a highly geometry dependent proportionality factor.Our linear solver can directly handle such a linear sliding relation.In order to solve the full nonlinear Eq. (B3) an iterative method is chosen that consecutively updates τ b until convergence is reached (for more details see analogue viscosity update described in Fürst et al., 2011).

Fig. 1 .
Fig. 1.Surface velocity magnitude for the IS standard interglacial equilibrium with continental outline in grey.The underlying geometry is obtained with the DR SIA model.(a) and (c) show the velocity magnitude for DR SIA while their higher-order equivalent from SR HO is presented in (b) and (d) for the same geometry.The lower two panels (c) and (d) show a close-up around Kangerlussauq.

Fig. 2 .
Fig. 2. Ice mass loss of different model versions on 5 km resolution for the three forcing scenarios: (a) MarAsl2, (c) TotAsl2, (d) MarCut started from IS while in panel (b) MarAsl2 is initiated from PS. Mass loss is determined by subtracting the background signal from a control experiment.The mass loss differs by no more than 20 % amongst the model versions and results are robust under different initial conditions (a) and (b).Each close-up zooms in on the last 5 yr of the experiments.The dashed blue (ME SIA) and green (SR HO) line represent the respective model version using a full nonlinear sliding law (cf.Appendix B).

Fig. 3 .
Fig. 3.Additional mass loss from changes in surface mass balance for the MarAsl2 experiment initiated from IS with respect to a control run.This mass loss component increases almost linearly and higher-order dynamics reduce the magnitude.

Fig. 4 .
Fig. 4. Initial response in surface velocities to the MarAsl2 perturbation for the three main model versions.Initiated from the 5 km IS, the geometry is still identical although the A sl forcing already operates.Shown are accelerations with respect to a control experiment run for each model version.In the DR SIA setup (a), the marginal accelerations are inherently confined to the forcing area.The close-ups of Jakobshavn Isbrae (d) and Kangerlussauq (e) allow the clear identification of the grid cells that experience A sl amplification.In the ME SIA (b, f, g) and SR HO setup (c, h, i), direct horizontal coupling is included and effects from both longitudinal stresses and transverse horizontal shearing are distinguishable.

Fig. 5 .
Fig. 5. Ice thickness evolution for the 5 km MarAsl2 experiment for all five model versions started from the IS state.The upper panels show the DR SIA ice thickness differences with respect to its control run after respectively 1, 25, 50 and 100 yr.Ice thickness changes in the other model versions are shown with respect to the dominant DR SIA differences.

Fig. 6 .
Fig. 6.Differences in ice velocity changes for the 5 km MarAsl2 experiment after 50 yr for ME SIA (a, c) and SR HO (b, d).Surface (upper panels) and basal velocity changes (lower panels) are shown with respect to the changes in the DR SIA velocity evolution.

Fig. 7 .
Fig. 7. Decomposition of changes in ice discharge for MarAsl2 into contributions from different evolution in (a) ice thickness alone, (b) velocity magnitudes alone and (c) their combination (cf.Eq. 1).Together with changes in surface mass balance all these curves add up to the total volume loss in Fig. 2a.

Fig. 8 .
Fig. 8. Spatially averaged velocity response in MarAsl2 for the marginal sliding region decomposed into its basal (a) and mean deformational component (b).Background velocities from an unforced control experiment are subtracted.Due to changes in the extent of the sliding area, both velocity components were smoothed using a 5 yr running mean.

Fig. 9 .
Fig. 9. Mean reaction time retrieved for each experiment: (a) MarAsl2 and (b) MarCut.The dot sizes indicate the resolution on which the experiment was conducted.The largest and smallest dot size represents the 5 and 20 km grid, respectively.Reaction times are calculated based on Eq. (2).This timescale reflects the immediate response to a stepwise perturbation rather than a full ice sheet equilibration of many thousands of years.

Table 2 .
Overview of model versions with different approaches to describe ice deformation and basal resistance.Note that the ME HO combination is omitted in this work.Details to each approach are described in Sect.2.3 and in Appendix B.

Table 3 .
Overview of the three perturbation experiments.