Modelling the dynamic response of Jakobshavn Isbræ , West Greenland , to calving rate perturbations

Jakobshavn Isbrae, one of Greenland’s major outlet glaciers, displayed rapid changes since the mid-1990s. Its floating ice tongue broke up around 1997, followed by a rapid calving front retreat over tens of kilometres, possibly linked to a prior warming of the ocean waters adjacent to the fjord. Parallel to this major retreat, a quasi-simultaneous process of acceleration and thinning of the glacier has been observed, currently making it a major contributor to eustatic sea-level rise. However, the causal interplay between the various factors involved has not yet been fully understood. Numerical studies of Jakobshavn Isbrae so far are either 2-D plan view ice flow models with a fixed calving front or 2-D flowline models, which have to parameterize lateral stresses. Hence the interaction between changes in calving front position and ice dynamics could not be studied consistently. To overcome this limitation, we implemented an implicit boundary tracking method in the Ice Sheet System Model (ISSM). This tool allows us to freely evolve the calving front by prescribing a calving rate, the ice front velocity being therefore the sum of ice velocity and calving rate. A suite of sensitivity experiments perturbing an initial steady-state calving rate has been performed to study its impact on the dynamics of the glacier. Our numerical results suggest a high sensitivity of the glacier dynamics to the applied calving rate. Changes in calving rate quickly affect upstream areas of the ice stream through a combination of changes in calving front position, ice velocity, thickness, grounding and ungrounding, and surface gradient change. Consequently, acceleration triggered at the calving front quickly affects the entire drainage basin. Moreover, the model results suggest that the ice stream does not recover from a short duration (about a year) calving rate perturbation over timescales on the order of a century. We present selected results of the sensitivity experiments to support the discussion clarifying the causes of the current changes occurring at this dynamic ice stream.


Introduction
Calving of icebergs is a major mean of ice discharge for marine terminating glaciers around the world.It accounts for about 50 % of the ice discharge of the Greenland and Antarctic Ice Sheets (Cuffey and Paterson, 2010;Rignot et al., 2013).Calving causes Introduction

Conclusions References
Tables Figures

Back Close
Full ice front retreat, which leads to reduced basal and lateral stress and results in upstream flow acceleration.Understanding calving dynamics remains challenging because of the diversity of factors that can trigger calving events.Bathymetry, tides and storm swell, as well as sea ice cover and temperatures of both sea water and air are possible factors influencing calving rates.However, their effect, their respective share and their interplay seem to vary from glacier to glacier, and is not well understood.Therefore, no universal calving rate parametrisation, which can be used in macroscopic ice sheet models, exists to date (Benn et al., 2007;Cuffey and Paterson, 2010).
In order to assess the impact of calving on the dynamics of an outlet glacier, we need to implement boundary conditions at the dynamically evolving ice front.This poses a second technical challenge as the ice front position needs to be tracked, and the related boundary conditions need to be updated accordingly.Addressing these issues is rather straightforward for 1-D-flowline or 2-D-flowband models (Nick et al., 2009;Vieli and Nick, 2011), where the ice front can be tracked along the floor line.However, this type of models lacks the consistent representation of mechanisms like lateral momentum transport and ice influx of converging ice streams, which have to be parametrised instead.This parametrisation may neglect feedback effects important for simulations on decadal to centennial time scales, e.g.catchment area entrainment (Larour et al., 2012a).
Therefore, it is desirable to implement a front tracking scheme in 2-D-horizontal models, which has been addressed by a few ice sheet models only, e.g.Winkelmann et al. (2011).Various approaches to model the evolution of the shape of ice exist.Explicit methods track the position of a set of points, which represent the ice front.They require a complex technical framework to allow for geometric operations like folding and intersection of the continuum boundary, tracking singularities in curvature, and determining the position of a point in space relative to the modelled continuum.On the other hand, a suite of implicit boundary tracking methods exists, e.g. the Level-Set Method (LSM) by Osher and Sethian (1988) implicitly by a contour, or "level-set", of a scalar valued "Level-Set Function" (LSF).It easily handles topological changes of the modelled continuum, like splitting and merging of its parts.It is well suited for the use in a parallel architecture, since it is based on a partial differential equation.This allows for application to large scale ice sheet simulations.The method is well established in Continuum Fluid Mechanics (Chang et al., 1996;Groß et al., 2006), and can be used with various numerical schemes, like the Finite Differences Method or the Finite Element Method (FEM).A LSM has been applied to ice modelling in test cases (Pralong and Funk, 2004), but not to real ice sheets yet.Jakobshavn Isbrae is a major marine terminating glacier in West Greenland, which drains about 6.5 % of the Greenland Ice Sheet (Zwally et al., 2011).It is characterised by two tributaries, which today terminate into an ice-choked fjord of about 30 km length (Figs. 1 and 2).The southern tributary exhibits high flow velocities, which are confined to a narrow, deep trough of about 5 km width.The trough retrogradely slopes inland to a maximum depth of about 1700 m below sea level (Gogineni et al., 2014).It discharges most of the ice of the drainage basin (Rignot and Mouginot, 2012).Observations have shown that the fast flowing areas of Jakobshavn Isbrae exhibit a weak bed with a considerable layer of temperate, soft ice (Lüthi et al., 2002).Most of these areas' horizontal motion is due to basal sliding and shear in this layer.A large fraction of the ice stream's momentum has to be transferred to the adjacent ice sheet by lateral stress.It is thus well justified to use the two-dimensional shelfy-stream approximation (SSA, MacAyeal, 1989) for model simulations of this glacier.
The basal topography of Jakobshavn Isbrae and its velocity field permit to roughly partition the catchment area into the confined, deep and fast-flowing troughs ("ice stream") and the surrounding parts of low velocity flow ("ice sheet").Those areas are separated by pronounced shear margins to either side of the ice stream.
Until the late 1990s, the glacier had a substantial floating ice tongue, which extended well into the ice fjord, and was fed by both tributaries.The position of the ice front remained fairly constant over time since the 1960s, and the glacier exhibited negligible seasonal variations in flow speed (Echelmeyer and Harrison, 1990).The glacier started Figures

Back Close
Full an ongoing phase of acceleration, thinning and retreat, that followed the breakup of its ice tongue.Seasonal variations in ice front position and flow velocity increased sharply (Joughin et al., 2004(Joughin et al., , 2008)).The glacier is currently far from equilibrium, making it one of the fastest ice streams in the world and a major contributor to global sea level rise (Howat et al., 2011;Joughin et al., 2014).Observations confirm that the current ice stream dynamics are mainly controlled by the position of the ice front (Podrasky et al., 2012;Rosenau et al., 2013;Moon et al., 2014).
Various hypotheses have been proposed to explain the mechanisms behind the change.All identify the breakup of the floating ice tongue as the initial triggering mechanism of this dramatic chain of events, but different possible mechanisms have been proposed to explain the sustained acceleration, thinning and retreat of the glacier.Studies by Joughin et al. (2012) and Habermann et al. (2013) propose loss of buttressing and changes in basal conditions as the main cause behind the ongoing acceleration.On the other hand, van der Veen et al. (2011) argue that the acceleration has to be accompanied by significant weakening of the lateral shear margins.Various modelling studies of the glacier, which use 1-D-flowline and 2-D-flowband models, project unstable retreat of the glacier along its southern trough for up to 60 km inland within the next century (Vieli and Nick, 2011;Joughin et al., 2012;Nick et al., 2013).Other modelling studies argue that this type of ice stream displays stable behaviour as long as it is being fed by the surrounding ice sheet (Truffer and Echelmeyer, 2003).However, numerical modelling efforts of Jakobshavn Isbrae in the horizontal plane so far lacked the representation of a dynamically evolving ice front.Hence, the hypotheses could not be tested in a satisfactory manner.
We present here a theoretical and technical framework for a LSM used to model the dynamic evolution of an ice front.This method tracks the ice front related boundary conditions and is a step towards a more physical representation of ice front dynamics in 2-D and 3-D ice sheet models.We implemented the method into the Ice Sheet System Model (ISSM, Larour et al., 2012b) perturbations in calving rate and conclude on its sensitivity to ice front retreat over the next 120 years.

Ice flow model
We employ the SSA on both floating and grounded ice.It neglects all vertical shearing but includes membrane stresses.The ice viscosity µ follows Glen's flow law (Glen, 1958): (1) Here, n = 3 is Glen's flow law coefficient, B the ice viscosity parameter, and εe the effective strain rate.We apply a Neumann stress boundary condition at the ice-air and ice-water interface, corresponding to zero air pressure and hydrostatic water pressure, respectively.A linear friction law links basal shear stress σ b to basal sliding velocity v b on grounded ice: where α denotes the basal friction parameter.We calculate effective basal pressure N assuming that sea water pressure applies everywhere at the glacier base.The ice thickness H evolves according to the mass transport equation: Here, v is the depth-averaged horizontal ice velocity, and a s and a b are surface and basal mass balance, respectively.We determine the grounding line position using hydrostatic equilibrium, and treat it with a sub-element parametrisation (Seroussi et al., Introduction Conclusions References Tables Figures

Back Close
Full 2014).We refer the reader to Greve and Blatter (2009), Cuffey and Paterson (2010) and Larour et al. (2012b), respectively, for further details on ice flow modelling and its implementation in ISSM.

Level-Set Method
Let Ω be a computational domain in two-or three-dimensional space, and ϕ a real, differentiable function on Ω × [0, ∞).We use ϕ to partition Ω at time t into three disjoint subdomains Ω i (t), its complement Ω c (t), and their common boundary Γ(t).We omit the time dependence of these sets in the following, unless stated otherwise.Let x be a point in Ω.If ϕ(x, t) < 0, then x belongs to Ω i .If ϕ(x, t) = 0, then x lies on Γ.If ϕ(x, t) > 0, then x belongs to Ω c .By construction, Γ, the 0-contour or "0-level-set" of ϕ, separates Ω i and Ω c .A particle at the boundary Γ moves with the boundary speed w .This motivates the "Level-Set Equation" (LSE): This Hamilton-Jacobi type partial differential equation describes the transport of the ice boundary across Ω.We need to provide an initial field ϕ 0 (x) = ϕ(x, t = 0) on Ω to solve Eq. (4).We can propagate the unit surface normal n on Γ onto Ω using the LSF by By definition, n always points into Ω c .For more details on the Level-Set Method and its applications, we refer to Osher and Sethian (1988) and Sethian (2001).
The boundary of an ice sheet evolves with a velocity, which is the sum of the ice velocity and an ablation velocity a = −a ⊥ n.The ablation rate a ⊥ = −a•n is non-negative Introduction

Conclusions References
Tables Figures

Back Close
Full Note that no limitations have been made with respect to the dimension of the problem in this section.Accordingly, the method could be applied to model the evolution of the glacier thickness and lateral extent simultaneously (Pralong and Funk, 2004).However, here we use the LSM to model only the horizontal extent of the ice sheet.Its vertical extent is described by the mass transport equation (Eq.3).We assume in the remainder of the article that lateral ablation occurs in the form of calving with a calving velocity c = −c ⊥ n, and that calving itself is a quasi-continuous process, consisting of frequent, but small calving events.With Eqs. ( 5) and ( 6), Eq. ( 4) becomes which is also known as "Kinematic Calving Front Condition" (KCFC, Greve and Blatter, 2009).The calving rate needs to be provided as a 2-D field on Ω, that can also vary with time.An example of a calving rate field can be seen in Fig. 4. The ice front Γ is advected at any time with the sum of the horizontal ice velocity and local calving rate values at the current ice front position (Fig. 3).We define the "calving flux" to be the ice volume calving at the ice front per unit time: Calving rate and velocity need to be given on the entire computational domain Ω in order to solve the KCFC.Since field variables like velocity and thickness are only Introduction

Conclusions References
Tables Figures

Back Close
Full available on Ω i , we need to propagate them continuously onto Ω c .We extrapolate such a field variable S onto Ω c by solving while keeping it fixed on Ω i .

Implementation
ISSM relies on the FEM to solve partial differential equations.It applies a Continuous Galerkin FEM using triangular (2-D) and prismatic (3-D) Lagrange finite elements, and uses anisotropic mesh refinement to limit the number of degrees of freedom while maximizing spatial resolution in regions of interest.
We discretise the KCFC (Eq.7) and extrapolation equation (Eq.9) using linear finite elements on the same mesh as the one used to model the ice dynamics.We stabilise both equations with artificial diffusion (Donea and Huerta, 2003), which after thorough testing proved to be the most robust stabilisation scheme.We integrate over time using a semi-implicit finite difference scheme.We solve the KCFC and the field equations for ice modelling in a decoupled fashion.The KCFC is solved first with input data from the last time step.We then update the numerical ice domain Ω h i using the new LSF, and update boundary conditions accordingly.We finally solve the momentum balance and the mass transport equation on Ω i .
The 0-level-set of ϕ does in general not coincide with the finite element mesh edges due to its implicit representation.The ice front Γ intersects a number of elements ("front elements") with a hyperplane, which divides them into an ice-filled and an ice-free domain (Fig. 3).Ice velocity and ice thickness exhibit a discontinuity at subgrid scale here, which we cannot resolve using the Continuous Galerkin FEM.This has various implications on the numerical level.When assembling the system stiffness matrices for ice modelling, exclusive integration over the ice-filled part of the element would be required.The stress boundary condition at the ice front would have to be applied at the Introduction

Conclusions References
Tables Figures

Back Close
Full intersecting hyperplane.Currently, ISSM is not capable of resolving those submesh scale processes.Therefore, we opt to either fully activate or deactivate an element at every time step.Only active elements are considered for the numerical discretisation of the respective field equations.We activate a mesh element if at least one of its vertices is in Ω i or Γ.
Then we consider the entire element to consist of ice.We flag the element as ice free if it lies entirely inside Ω c , and it is deactivated at this time step.As a consequence, the numerical ice front Γ h runs along mesh edges, and updates of Γ h occur in a discontinuous manner.We apply the stress boundary condition along Γ h for numerical consistency.The strongly viscous ice rheology efficiently dampens effects of these discrete boundary condition updates.Ice front normals on Γ and Γ h may differ significantly in direction.However, stress components tangential to n cancel out along Γ h , so that the integrated stress exerted at the ice front is close to the one applied along Γ.For all further calculations where a normal is involved, like extrapolation, the normal to the LSF (Eq.5) is used.
The ice domain Ω i is by definition a subset of the numerical ice domain thus displaced further downstream than Γ.This may lead to slightly higher resistive lateral stress at the ice front, whose magnitude depends on the excess ice area of the intersected front element and the front geometry.We choose a fine mesh resolution in the vicinity of the ice front to limit this effect.
We extrapolate the ice front thickness onto the ice-free vertices of the front elements using Eq. ( 9).This yields realistic ice thickness and ice thickness gradients across the front elements, that would otherwise lead to overestimated driving stress and underestimated water pressure at the ice-ocean interface.If not corrected, those two effects unrealistically increase ice velocities at the ice front, which then feed back into the mass transport and LSM schemes.
We present an idealized test-setup for validation of the LSM-module in Appendix A. It shows that the ice margin is advected with correct speed.A small error is due to the underlying unstructured mesh used, but cancels out over time.The module re-Introduction

Conclusions References
Tables Figures

Back Close
Full quires additional computational effort for the extrapolation of field variables, to solve the KCFC, and for extra iterations of the momentum balance solver, since the stress boundary condition at the ice front does often change.Its amount depends on the flow approximation used and especially on whether the model setup is close to a stable configuration or not.Using the SSA approximation, the additional computational cost amounts to up to 25 %, of which 11 % is caused by the solution of the LSM-module for this experiment.
3 Data and model setup

Jakobshavn Isbrae model setup
We use Jakobshavn Isbrae's drainage basin from Zwally et al. (2011) to generate a 2-D-horizontal finite element mesh with element size varying from 500 m in the fjord and areas of fast flow to 10 km inland (Fig. 4).We choose this high mesh resolution to minimise ice front discretisation errors, and to resolve the fjord and the deep trough accurately in the model.The resulting mesh has about 10 000 vertices and 19 000 elements.Due to high flow velocities, the numerical Courant-Friedrich-Lewy condition (Courant et al., 1928) dictates a time step length on the order of days.We use ice bedrock topography from Morlighem et al. (2014), derived using a mass conservation approach (Morlighem et al., 2011).The ice surface elevation is taken from GIMP (Howat et al., 2014).Ice thickness is the difference between ice surface and ice base elevation.Bathymetry of the ice-choked fjord of Jakobshavn Isbrae is difficult to measure and currently poorly known.As a first order estimate, we apply a parabolic profile of 800 m depth along the ice fjord, fitted via spline interpolation to known topography data.We rely on Ettema et al. (2009)  of the model simulations.Basal mass balance is set to zero and no thermal model is run.All these forcings are kept constant over time.
We infer a basal friction coefficient α in Eq. ( 2) using an adjoint-based inversion (MacAyeal, 1993;Morlighem et al., 2010) of existing InSAR-derived surface velocities from 2009 (Rignot and Mouginot, 2012).In regions like the ice fjord, where there is no ice today, we apply an averaged value of α = 30 a 1/2 m −1/2 .At the margins of the computational domain we prescribe zero horizontal ice velocities in order to prevent mass flux across this boundary.α is kept fixed over time for all model simulations.Inconsistencies in model input data cause sharp readjustments of the glacier state at the beginning of each simulation, which would make it difficult to distinguish between such effects and those of the applied forcing (Seroussi et al., 2011).Therefore, we relax the ice surface prior to the experiments using a fixed ice front, whose position and orientation is set to the mean annual ice front position of 2009.Since the glacier in this configuration is far from steady state, model relaxation causes considerable thinning across the glacier's catchment area.In order not to deviate too much from present day's geometric setting we choose a 100 year relaxation time interval.Note that the grounding line retreats during the relaxation due to dynamic thinning, so that the glacier forms a new floating ice tongue.This ice tongue extends about 15 km to a topographic high in the southern trough and 3 km into the northern one (Fig. 4).The relaxed geometry constitutes the initial state for our experiments.

Description of experiments
For simplicity, frontal ablation occurs in the experiments exclusively in the form of calving.We let relaxation run, continued onto Ω c (Fig. 4).l is a continuous function, which is 1 in areas where the bedrock lies below −300 m, and linearly drops to 0 in areas of positive bedrock elevation.It prevents calving onto areas with a glacier bed above sea level, motivated by observations from tidewater glaciers (Brown et al., 1982).We expect this calving rate to yield a stationary ice front, if applied to a geometry in steady state.We scale c ⊥ 0 over time with a scaling function s = s(t), which allows for the representation of seasonal calving cycles, and a perturbation function p = p(t) to modify the calving rate for some period of time, respectively.The applied calving rate then is c We perform three suites of experiments in order to analyse the impact of the calving rate on the glacier's dynamics.The ice front is now allowed to freely evolve in response to c ⊥ .All experiments run for a total time of 120 years.
In experiment A, we keep calving constant over time, i.e. we let both . This experiment, although not physically motivated, is used to evaluate whether a stable ice front position can be reached using the LSM.
In experiment suites B and C, we mimic seasonal calving by scaling c ⊥ 0 over time by a factor s(t) = max(0, π sin(2π (t/L − φ 0 ))), with phase shift φ 0 = 4/12 and period L = 1 a.We perturb the calving cycle during a limited perturbation duration ∆t with a perturbation strength p 0 ≥ 0: can be compared to.We introduce P = (1 − p 0 ) • ∆t, which is a measure of the timeintegrated deviation in applied calving rate relative to B1.

Results
We present ice front positions for several experiments in Fig. 5.Under constant calving rate forcing, the ice front essentially remains at a stable position after minor readjustments in the first decade of the simulation.In experiment A, the ice front undergoes gradual retreat over time.When we perturb the calving rate, the ice front migrates.
Increased calving rates cause ice front retreat: higher calving rates lead to larger retreats.The retreat is highest in areas of fast ice flow, but strongly decreases towards the ice stream margins.This yields the characteristic concave shape of a retreating ice front.Resulting ice fronts positions and their characteristics do agree well with ice fronts obtained from observations (Fig. 2).Ice front positions of all experiments B and C readvance after the perturbation interval to a position which is indistinguishable from the one of B1 in one to two decades.
Figure 6 shows ice velocity, geometry and strain rates along two tracks, which go along and across the southern trough, respectively (Fig. 4).We observe that the ice front position, its ice velocity and thickness as well as the position of the grounding line form an ensemble, which quickly adjust to each other under constant mean annual calving rates to reach a characteristic "ice front configuration".The variables of the ice front configuration mirror the minor adjustment of the ice front position prior to the perturbation.Only the ice thickness in the floating part decreases up to 100 m once the ice front is allowed to freely adjust.Enhanced calving causes the modelled ice stream to rapidly thin, accelerate and retreat.However, the ice front can be very thick at the end of the calving season, as the glacier has retreated back into areas of high ice thickness.
Ice velocities respond to this change by increasing over time as the ice front retreats and peaking at the point of furthest retreat.The response is stronger with larger and longer perturbations, as expected.The acceleration of the ice stream extends well up-Introduction

Conclusions References
Tables Figures

Back Close
Full stream to areas of grounded ice.Both thinning and acceleration mainly occur in the ice stream itself, and diffuse out to the surrounding ice sheet in a strongly dampened fashion.These thinning and acceleration patterns increase surface and velocity gradients along and across the ice stream, especially in the shear margins.Here, the ice stream acceleration gradually increases effective strain rates up to a factor of 4 in experiment C4, which corresponds to a drop in viscosity of 60 % (Eq.1).This weakens the mechanical coupling between the ice stream and the surrounding ice sheet substantially.The grounding line position lies downstream of a local high of the basal topography at any given time of the experiment, even if the general slope of the basal topography is retrograde.It is hardly affected by small calving rate perturbations.However, large perturbations cause discontinuous migration of its position to the next topographic high upstream.Their discrete locations form the set of stable grounding line positions along the deep part of the tributaries.Discontinuous retreat of the grounding line causes drastic, but short-lived flow acceleration, which indicates that the basal topographic highs act as pinning points, even if basal friction coefficient α and effective basal pressure N in the troughs are low.Timelines of ice front position, ice velocity and thickness at the ice front, as well as effective strain rate and position of the first grounded point along-trough are shown in Fig. 7.They allow us to study the interannual behaviour of the ice front configuration in more detail.All shown variables reflect the characteristics of the applied calving rate forcing.They stay constant in experiment A, and oscillate periodically with a constant amplitude around an annual mean value in experiment suites B and C.
Prior to the perturbation, the ice front position varies by ±3 km along-stream.Ice velocities and front thickness act in phase to the ice front position, while the response of strain rates at the grounding line is slightly delayed.Ice velocities vary by ±20 %, which corresponds to about ±2 km a −1 , ice thickness by ±13 %, or ±100 m, and effective grounding line strain rates by ±7 %, or ±0.1a −1 .The grounding line position is stable at kilometre 29 of the along-trough profile.In response to a two-fold increase in calving rate, the ice front retreats initially at an average rate of 4.5 km a −1 (Fig. 7d).Introduction

Conclusions References
Tables Figures

Back Close
Full This retreat rate decreases to zero for longer perturbations, so that the ice front approaches a new stable position upstream.Variation of the ice front position during the retreat doubles to up to ±6.5 km.Average ice velocity increases by about 10 %, but its variation also almost doubles to ±38 %.Occasional velocity spikes occur related to ungrounding from pinning points.The ice front thins on average by a third, with oscillations up to ±75 %.This high thickness variation is due to the calving back into areas of thick ice in summer followed by stretching and thinning during ice front advance in the late winter season.Variation of the grounding line position remains small during the perturbation, disrupted only by its discontinuous retreat for large perturbations.
For small perturbations (Fig. 7c), variation of effective strain rates here quadruples to ±25 %.Once the calving rate perturbation ceases, the ice front configuration displays a striking reversibility in all its variables.Finally, Fig. 8 shows the evolution of the total ice volume as modelled against control run B1.The mean annual volume change due to ongoing geometry relaxation of experiment B1 is −22.8 km 3 a −1 , varying by ±10.5 km 3 per year.Experiment A undercuts this value by −0.4 km 3 a −1 on average related to the gradual retreat of its ice front.Enhanced calving causes additional volume loss proportional to perturbation measure P .It amounts to −35.7 km 3 a −1 in the first year of perturbation of experiment C8, but decreases with time, as the ice front thins and retreats into areas of lower calving rates.
Those numbers agree well with recent ice discharge observations (Howat et al., 2011).advances during the perturbation, creating a convex ice tongue.Meanwhile, the ice stream decelerates and thickens, causing grounding line advance and volume increase.After the perturbation period, the glacier calves back into a state slightly thicker and faster than the one of B1.

Discussion
The characteristics of the applied calving rate determine the behaviour of the ice front.
The ice front will reach a stable configuration if the calving rate exceeds the ice velocity from some point along-stream onward.It will either retreat back or advance to this point.If no such point exists, the ice front will undergo either unconfined advance or retreat, depending on the magnitude of the calving rate.
Adjustments of the ice front configuration to the current calving rate change the local driving stress, which lead to corresponding changes in thickness and surface gradient, and vice versa.The mass transport mechanism allows the thickness change to diffuse upstream.In case of ice front thinning and surface steepening, the increased driving stress leads to enhanced mass flux, causing further upstream thinning, steepening and acceleration.
We observe multiple secondary feedback mechanisms in addition to the mass transport mechanism, that determine the way the ice sheet adjusts to calving rate forcings.Thinning of the ice stream leads to reduced basal effective pressure and grounding line retreat, which both reduce basal drag significantly.Detachment of the base from pinning points leads to short-lived, but drastic increases in ice flux.These mechanisms are qualitatively the same as the ones described in Vieli and Nick (2011) and Joughin et al. (2012).
The ice front lengthens during its retreat, which increases the total driving stress exerted by the ice front.Side arms of the main tributaries are now calving directly into the fjord.This increases calving flux Q cf , which additionally thins and steepens areas in the immediate ice front vicinity.

Conclusions References
Tables Figures

Back Close
Full The steepening of the surface slopes across the shear margins caused by ice stream thinning results in increased driving stress, which increases lateral ice inflow into the ice stream.This in turn limits surface steepening and grounding line retreat along the trough.The feedbacks described above cause a net rise in volume flux towards the ice front and enable quick adjustment of the glacier in response to changes of the ice front configuration.
Acceleration of the ice stream and surface steepening in its vicinity strongly increase strain rates at the shear margins and the grounding line.The non-linear rheology of ice softens these areas, allowing the ice stream to accelerate further, and to soften its shear margins more and more.This positive feedback is only controlled by the rate at which ice is able to enter the ice stream to sustain the driving stress.We consider this mechanism to be substantial for sustaining the longer term ice stream acceleration, since a large fraction of the ice stream's driving stress is balanced by lateral stress.This corroborates force balance arguments produced earlier by van der Veen et al. (2011).
Reduced lateral momentum transfer decreases lateral mass influx.This makes the surrounding ice sheet less sensitive to short periods of enhanced calving.For this reason, calving induced thinning mainly occurs in the ice stream, which limits the rate at which additional ice can be discharged.
Conversely, in case of no seasonal calving cycle of experiment A, mechanical coupling between the ice stream and ice sheet is higher than compared to experiment B1.Ice stream velocity is therefore lower, which causes net ice front retreat and corresponding volume loss.This illustrates that ice sheet models which do not include a seasonal calving cycle may overestimate mass loss of glaciers.Moreover, we do advise against the use of ice models, which do not incorporate a dynamically evolving ice front nor the related lateral effects, for projections of future contributions to global sea level rise on decadal to centennial time scales.
Response mechanisms not covered here will likely include feedbacks in damage mechanics and thermodynamics due to the increased strain rates.During longer per-Introduction

Conclusions References
Tables Figures

Back Close
Full turbations, ice surface lowering will probably affect the surface mass balance and the drainage basin outline.The reversibility of the ice front configuration after the calving rate perturbation is a robust feature across all experiments.We see the short duration of the perturbation, the prescribed calving rates and the geometric setting of the ice stream to be responsible for this behaviour.Volume change relative to experiment B1 is mainly due to the migration of the ice front and the according change in lateral ice sheet extent.However, the additional ice discharge never exceeds 0.1 % of the total glacier volume in the experiments shown here.This volume loss can be easily balanced by the vast surrounding ice sheet.Since calving rates are not coupled to ice dynamics, the ice front configuration is able to quickly recover once this forcing is set back to its initial value.
The modelled glacier response due to enhanced calving is in good agreement with observations.However, the reversibility of the ice front position is in stark contrast to Jakobshavn Isbrae's ongoing observed retreat.Therefore, the actual calving rates must have stayed increased, as our results suggest that the glacier would have readvanced otherwise.

Conclusions
We presented the theoretical framework for coupling a LSM to ice dynamics and implemented it into the ISSM.The LSM proved to be a robust method for modelling the dynamic evolution of an ice front subject to ablation.We applied the technique to Jakobshavn Isbrae using prescribed calving rates.We find its dynamics to be highly sensitive to the applied calving rate, which agrees well with observations.
The calving rate perturbation strongly affects the ice flow inland through various mechanisms.Changes in calving rate cause ice front migration and alter its ice discharge.The resulting thickness change at the ice front diffuses upstream through a coupling between the stress balance and mass transport mechanism.Ice stream thinning reduces basal drag by means of grounding line retreat and reduction of effective basal pressure.Rheological shear margin weakening caused by ice stream acceleration decreases lateral drag resistive to flow.These two positive feedback mechanisms are able to sustain prolonged acceleration of the ice stream, and are only controlled by the rate at which ice can enter the ice stream.However, the vast surrounding ice sheet is barely affected by short periods of enhanced calving.It stabilises the ice stream and allows for quick reversibility of the ice front configuration through lateral mass influx and stress transfer once we set the calving rates back to their initial value.The importance of the ice front position and lateral effects for the behaviour of this type of glacier lets us advise against the use of models, which do not represent these mechanisms, for projections of future ice discharge.
The method presented here is able to close a major gap present in ISSM and various other ice sheet models.It enables a multitude of future studies in combination with e.g.thermodynamics and damage mechanics, and will improve our understanding of outlet glaciers.

Appendix: Test-setup
We present a simple test setup to validate the LSM.Let Ω be a 50 km square with an initial LSF as where x 0 = (25, 25) km and R = 12.5 km, so that our initial 0-level-set describes a circle in the middle of the domain.We prescribe a constant velocity of v = 1 km a −1 • (cos(π/4), sin(π/4)) everywhere.We advect ϕ 0 over 10 years with time steps of 0.1 a, and keep track of its 0-level-set.
We plot the 0-level-set position at the beginning every year in Fig. A1.The initial circular shape is preserved and advected with constant speed.We measure the advection speed of the 0-level-set along the diagonal marked in white in Fig. A1. Figure A2 shows the standard deviation of its numerical error relative to the prescribed advection speed Introduction

Conclusions References
Tables Figures

Back Close
Full for different element sizes.The numerical error is mainly related to the interpolated representation of the curved shape on an irregular mesh.This causes the front velocity to vary around the mean value, which indeed equals the prescribed velocity.The standard deviation of the error linearly decreases with mesh width, and drops below 1 % for elements sizes below 0.5 km.We therefore recommend using an element size of below 1 km in the ice front vicinity for ice sheet simulations.Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | The LSM represents the continuum boundary Figures Back Close Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | , a parallel, state-of-the-art FEM ice sheet model, and apply it here to Jakobshavn Isbrae in order to model the dynamic response to Figures Back Close Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | in compliance with current glaciology terminology, and is the difference between ice boundary velocity w and ice velocity v projected along n a ⊥ = (v − w ) • n. (6) It follows that the ice boundary is stationary (w • n = 0) if and only if a ⊥ = v • n, i.e. the ablation rate matches the ice velocity perpendicular to the ice boundary.The ablation rate a ⊥ can be an input taken from observations or a suitable parametrisation.
Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | to force the surface mass balance.Their surface temperatures are used to calculate the ice viscosity parameter B following an Arrhenius relationship (Cuffey and Paterson, 2010) once at the beginning Discussion Paper | Discussion Paper | Discussion Paper | 10) as a basic calving rate estimate, motivated by the small observed angle between v and n at the ice front.Here, v 0 denotes the velocity field at the end of the geometry Discussion Paper | Discussion Paper | Discussion Paper | perturbation at t 0 = 20 a for all experiments.In experiment suite B we choose ∆t = 1 a and vary p 0 in steps of 1 from 0 to 4. On the other hand, in experiment suite C we keep p 0 = 2 fixed, and choose ∆t as 2, 4 and 8 years.We use the notation B< p 0 > and C< ∆t > to identify the single experiments, e.g.B2 for experiment B with perturbation strength p 0 = 2. Experiment B1 represents in this experiment suite the case of unperturbed periodic calving.It is used as a control run the other experiments Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | After the perturbation, all modelled glaciers recover 40-60 % of the volume difference in the first decade, because during the readvance the thin ice front strongly reduces the calving flux.The calving flux converges to the value of experiment B1 once the ice front configuration approaches the state of the control run.The ice volume now has to be replenished by excess surface accumulation, which takes longer than the 120 years of simulation time presented here.Calving rate perturbations therefore leave a lasting imprint on the ice sheet.Results of experiments of decreased calving are not shown here.They exhibit the exact inverse pattern to the one described above for enhanced calving: the ice front Figures Back Close Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |