Greenland ice sheet contribution to sea-level rise from a new-generation ice-sheet model

. Over the last two decades, the Greenland ice sheet (GrIS) has been losing mass at an increasing rate, enhancing its contribution to sea-level rise (SLR). The recent increases in ice loss appear to be due to changes in both the surface mass balance of the ice sheet and ice discharge (ice ﬂux to the ocean). Rapid ice ﬂow directly affects the discharge, but also alters ice-sheet geometry and so affects climate and surface mass balance. Present-day ice-sheet models only


Introduction
The currently observed acceleration of mass loss from the Greenland ice sheet (GrIS, Rignot et al., 2011;Schrama and Wouters, 2011;van den Broeke et al., 2009;Wouters et al., 2008) is a concern when considering its possible contribution to future sea-level rise (SLR).Approximately 60 % of the acceleration rate in mass loss from the GrIS has been attributed to a change in the surface mass balance (SMB, Rignot et al., 2011;van den Broeke et al., 2009).However, several studies have revealed a dynamic response of the ice sheet, in which acceleration and thinning of most outlet glaciers are shown to be responsible for a substantial increase in ice discharge (Howat et al., 2007;Pritchard et al., 2009;Joughin et al., 2010;Moon et al., 2012).These studies show a high spatial and temporal variability in glacier acceleration, suggesting that simple extrapolation of the recent observed trends cannot be justified, and realistic projections of the contribution of GrIS to SLR on decadal to century time scales must be derived from the forecasts of verified ice-flow models driven by the most reliable projections of climatic (atmosphere and ocean) forcing.
The flow of ice in an ice sheet is characterised by a very low Reynolds Number and is governed by the Stokes equations (e.g.Greve and Blatter, 2009).Outlet-glaciers dynamics are strongly controlled by basal and seaward boundary conditions.These boundary conditions have recently been altered by ongoing climate change: increased surface runoff filling the crevasses can result in a softening of the outlet glaciers lateral margins through cryo-hydrologic warming and hydraulic weakening of ice (Van der Veen et al., 2011) and can enhance basal lubrication by reaching the bed through moulins (Zwally et al., 2002); ocean warming and processes happening at the front have likely triggered the recent acceleration of numerous outlet glaciers by reducing the back-stress at the front as their floating tongues thin and/or retreat (Howat et al., 2007).
Despite the recent efforts to model these processes (Nick et al., 2009;Schoof, 2010), incorporating them and validating the results produced has not been the focus of most modellers running continental scale models (Vaughan and Arthern, 2007): most of those ice-sheet models were primarily designed to run over glacial cycles and consequently did not need to reproduce the decadal to annual sensitivity that recent observations have highlighted, and which are likely to be significant in decade-to-century projections (Truffer and Fahnestock, 2007).
The lack of skill of the ice-sheet models used for the Fourth Assessment of the Intergovernmental Panel on Climate Change, was one of the reasons behind the statement that a poor understanding of the importance of dynamic changes has limited our ability to put an upper bound on the contribution of ice sheets to sea level by 2100 (Solomon et al., 2007).Since then, the inherent fundamental limitations that prevent proper modelling of the ice discharge, have been identified: i.Most of the observed changes are located on narrow outlet glaciers and the dominant source of increased discharge is the combined contribution of many small glaciers (Howat et al., 2007;Moon et al., 2012) that cannot be captured individually by models typically running with grid resolutions from 5 km to 15 km; ii.Low order approximations of the Stokes equations do not hold for these outlet glaciers where the scale of horizontal variations in basal topography and friction are of the same order as the ice thickness (Pattyn et al., 2008;Morlighem et al., 2010); iii.Inverse methods, by offering a robust framework to combine information from the model and the data, are essential in modelling real systems but remain underused in glaciology.As a consequence, several parametrisations remain very poorly constrained limiting our ability to reproduce the current state of the ice sheet.The most uncertain parametrisation is the drag exerted on the ice by the underlying bed, this can vary by several orders of magnitude depending on the bedrock roughness and water pressure (Jay-Allemand et al., 2011).
We have developed a new generation continental scale icesheet model (Little et al., 2007;Alley and Joughin, 2012) that overcomes these difficulties: by employing parallel computing and the Elmer/Ice code (http://elmerice.elmerfem.org/),we solve the full system of Stokes equations over the entire GrIS (Sect.2.1).We use an anisotropic mesh-adaptation technique (Sect.2.2) to distribute the discretisation error equally through the entire domain (Frey and Alauzet, 2005;Morlighem et al., 2010).For the construction of the initial state, we use two inverse methods (Arthern and Gudmundsson, 2010;Morlighem et al., 2010) to constrain the basal friction coefficient field from observed present-day geometry and surface velocities (Sect.3.1).While other recent models have included some similar features (Price et al., 2011;Seddik et al., 2012;Larour et al., 2012), we present the first model to use all three developments simultaneously to produce prognostic simulations.Due to ice flux divergence anomalies (Seroussi et al., 2011) caused by the remaining uncertainties in the model initial conditions, the free surface is then relaxed for 50 yr (Sect.3.2).From this initial state, we run sensitivity experiments of the GrIS dynamical response over one century to perturbations in climate and basal lubrication (Sect.4).In conclusion, we discuss how these results can be interpreted in terms of possible bounds for the future ice-sheet mass loss and contribution to SLR by 2100.

Equations
We consider a gravity-driven flow of incompressible and non-linearly viscous ice flowing over a rigid bedrock.
The constitutive relation for ice is assumed to be a viscous isotropic power law, called Glen's flow law in glaciology (Glen, 1955): where τ is the deviatoric stress tensor, εij = (u i,j + u j,i )/2 are the components of the strain-rate tensor, and u is the velocity vector.The effective viscosity η is expressed as where εe = εij εij /2 is the strain-rate second invariant, E is an enhancement factor, A(T ) is the rate factor function of the temperature T relative to the temperature melting point following an Arrhenius Law: In Eq. ( 3), A o is the pre-exponential factor, Q is an activation energy, and R is the gas constant.
The ice flow is computed by solving the Stokes problem with non-linear rheology, coupled with the evolution of the upper free-surface, summarised by the following field equations and boundary conditions: On the domain , Eq. ( 4) expresses the conservation of mass and the conservation of momentum.The Cauchy stress tensor σ is defined as σ = τ − pI with p the isotropic pressure.The gravity vector is given by g = (0, 0, −g) and ρ i is the density of ice.Equation (5) expresses the evolution of the upper free surface s with a prescribed net accumulation rate a s (x, y, t).We note ∂ i z the partial derivative of upper surface z s (x, y, t) with respect to the horizontal dimension i = (x, y).A proper treatment of grounding line dynamics has been developed for three-dimensional full-Stokes simulations (Favier et al., 2011) but remains computationally challenging if applied to a whole ice sheet.Here, ice shelves and grounded ice are not treated differently and the lower surface elevation z b is fixed in time.
On the boundaries, n and t are the normal and tangential unit vectors.The upper surface s is a stress-free surface (Eq.6).A linear friction law is applied on the lower surface b (Eq.7).On the lateral boundary l , the normal component of the stress vector is equal to the hydrostatic water pressure exerted by the ocean, with ρ w the sea water density, where ice is below sea level l w , and is equal to zero elsewhere (Eq.8).The remaining components of the stress vector are null (Eq.8).
During transient simulations, the mesh is deformed elastically to follow the free-surface elevation with the constraint that the minimum thickness is equal to h min = 10 m.The mesh nodes are not allowed to move horizontally.As a consequence, for outlet glaciers, the calving front position is fixed in time, however, land terminating glaciers can retreat inland but not advance beyond the initial footprint (Sec.2.2).The

Parameters
Values Units ice-sheet volume is computed as the integral of the basis functions over the whole mesh.The ice discharge is computed as the ice flux through the lateral boundary l as As the glaciated area (i.e. the area where ice thickness is larger than h min ) can change with time, SMB is computed from Eq. ( 5) as Values of parameters prescribed in this study are presented in Table 1.A proper initialisation of the temperature field would require a spin-up over at least a full glacial cycle, which remains too computationally expensive with our model.In our application, the initial field is bi-linearly interpolated on the finite element mesh from the temperature field computed with the shallow ice model SICOPO-LIS (Greve, 1997) after a paleo-climatic spin-up as in Seddik et al. (2012).In order to avoid an initialisation shock, as these two models do not have the same physics, the ice temperature field is then kept constant.The low conductivity of ice means that this assumption remains valid for the short time scales considered here.Because the ice-sheet basal temperature conditions remain uncertain, basal sliding is not restricted to areas where the ice temperature has reached the pressure melting point and an optimal friction coefficient β in Eq. ( 7) is inferred by inverse methods at each node to reproduce the observed surface velocity field (Sect.3.1).
All the equations presented above are solved using the ice flow model Elmer/Ice, based on the finite element code Elmer (http://www.csc.fi/elmer).Details of the numerical methods used to solve Eqs. ( 4) and ( 6) and validation through established benchmark experiments can be found elsewhere (Gagliardini and Zwinger, 2008).We give a more detailed discussion in the implementation of the non-penetration condition, i.e the zero-normal velocity in Eq. ( 7), as the results on a rough topography will depend on the choice of the implementation strategy.Using finite elements, the normals are uniquely defined element-by-element using the nodes coordinates and elemental shape functions.At a given node, the discontinuity of the normal direction to the elements sharing the node can be large if the typical length scale of the basal roughness is smaller than the mesh size.
The non-penetration condition (Eq.7) can be enforced weakly by elements under the form of a Robin condition as where is an arbitrary large positive number.If the discontinuity of the normal direction, n, is large between adjacent elements, Eq. ( 11) will lead to u = 0 at the node shared by such elements.Enforcing Eq. ( 7) as a Dirichlet condition implies a linear combination of the three components of the velocity vector in the general case where n is not aligned with one of the coordinate axis.This can be done as in Leng et al. (2012) by using a local coordinate system with one direction aligned with the normal direction.A choice needs to be made on the definition of this normal direction as it is not uniquely defined, in general, at a given mesh node.The most immediate possibility is to take the average of the normal direction to the elements sharing the node.The mass flux through an element, computed as the integral of the scalar product of u with n, is then non null in general.This solution has been adopted in the following simulations.An alternative, newly implemented in the model, is to use mass consistent normals where the normal at a given node is a weighted average of the normal direction to the elements sharing the node.The weights are the nodal basis functions.This choice will insure that the zero-flux condition will be preserved globally, the non-null fluxes through the elements cancelling each other.However, it should be noted that this choice for the definition of the normal direction at the mesh nodes can be flawed with quadratic elements (Walkley et al., 2004).

Mesh construction
Anisotropic mesh adaptation is now widely used in numerical simulations especially with finite elements, as it allows to refine the mesh where needed to capture the flow features within a certain accuracy without increasing the computational cost excessively.The method is generally based on an estimation of the interpolation error used to adjust the mesh size so that the discretisation error is equally distributed over the whole domain.It can be shown that an estimate of the interpolation error induced by the meshing is obtained from the Hessian matrix of the modelled field, allowing to define an anisotropic metric tensor at each node (Frey and Alauzet, 2005).As the target surface velocity field is known, our metric tensor is constructed from the Hessian matrix of observed surface speed (Joughin et al., 2010) shown in Fig. 1a.Starting from a first regular mesh of the 2-D footprint of the GrIS based on Delaunay triangulation, we use the freely avail-  able anisotropic mesh adaptation software YAMS (Frey and Alauzet, 2005) to optimise the mesh sizes according to the given metric map.
The resulting mesh is depicted in Fig. 2. Mesh sizes decrease from 40 km in the central part of the ice sheet to 1 km in the outlet glaciers.The 2-D mesh is then vertically extruded using 16 layers.The resulting 3-D mesh is composed of 417 248 nodes and 748 575 wedge elements.
The bedrock and surface topography are taken from the freely available SeaRise 1 km present-day data set (http://tinyurl.com/srise-umt) and are based on the Bamber et al. ( 2001) digital-elevation models where new data have been added on three of the main outlets (Jakobshavn Isbrae, Helheim and Kangerdlugssuaq).There is no freely available product for current ice margin position and the initial 2Dfootprint has been constructed from a 0-m thickness contour.

Initial state
As ice-sheet responses include long time scales (multicentury), forecasting change on decadal-to-century time scales is essentially a short-term forecast.As such simula-tions are sensitive to the initial state, simulating the present conditions of the ice sheet is crucial (Vaughan and Arthern, 2007).Available observations of the current state of the ice sheet include the ice-sheet geometry (bedrock and freesurface elevations, e.g.Bamber et al., 2001), surface velocities (e.g.Joughin et al., 2010) and rate of change of the surface elevation (e.g.Pritchard et al., 2009).If time series of these observations are available for the last decade, observed changes in velocity and surface elevation are certainly the results of transient boundary forcing that are still not fully understood and modelled.Moreover, while the model developed by Heimbach and Bugnion (2009) offers the ability to do transient data assimilation, inverse methods applied to full-Stokes ice-sheet modelling are currently restricted to diagnostic simulations, hence limiting the ability to assimilate time series.In this application, we compare two recently developed inverse methods to constrain the basal friction coefficient field (β in Eq. 7) from a given geometry and surface velocity field, considered as representative of present-day conditions (Sect.3.1).Due to ice flux divergence anomalies caused by the remaining uncertainties in the model initial conditions, the free-surface-elevation rate-of-change F. Gillet-Chaulet et al.: GrIS contribution to sea-level rise computed in a diagnostic model is non-physical (Seroussi et al., 2011), and the available observations are not used here to constrain the model.The free surface is then allowed to relax compared to the observed surface for a period of 50 yr (Sect.3.2).The performance of the model to simulate present-day conditions is then assessed by comparison of the estimated ice discharge for the main outlet glaciers (Rignot and Kanagaratnam, 2006).

Inverse methods
Two variational inverse methods (Arthern and Gudmundsson, 2010;Morlighem et al., 2010) are used to infer the basal friction coefficient field β(x, y).Both methods are based on the minimisation of a cost function that measures the mismatch between modelled and observed velocities.The two methods are briefly outlined below.

Robin inverse method
The method, detailed in Arthern and Gudmundsson (2010), consists in solving alternatively the Neumann-type problem, defined by Eq. ( 4) and the boundary conditions in Eqs. ( 6) and ( 7), and the associated Dirichlet-type problem, defined by the same equations except that the Neumann uppersurface condition (Eq.6) is replaced by a Dirichlet condition where observed surface horizontal velocities are imposed.
The cost function that expresses the mismatch between the solutions of the two models is given by where superscripts N and D refer to the Neumann and Dirichlet problem solutions, respectively.The Gâteaux derivative d β J o of the cost function J o with respect to the basal friction coefficient β for a perturbation β is given by where the symbol |.| defines the norm of the velocity vector and b is the lower surface.
Note that this derivative is exact only for a linear rheology and thus Eq. ( 13) is only an approximation of the true derivative of the cost function when using Glen's flow law (Eq.1) with n > 1 in Eq. (2).

Control inverse method
The control method has been introduced by MacAyeal (1993) and recently applied to full-Stokes ice flow modelling by Morlighem et al. (2010).The method relies on the computation of the adjoint of the Stokes system.As in Morlighem et al. (2010), we assume that the stiffness matrix of the Stokes system is independent of the velocity and thus self adjoint, which is valid only for Newtonian rheology, i.e. when n = 1 in Eq. ( 2).
As the direction of the ice velocity is mainly governed by the ice-sheet topography, we disregard the error on the velocity direction and the cost function is expressed as the difference between the norm of the modelled and observed horizontal velocities as where u obs are the observed velocities and subscript H refers to the horizontal component of the velocity vector.
The Gâteaux derivative is obtained by where λ is the solution of the adjoint system of the Stokes equations.
Again, this derivative is exact only for a linear rheology and thus is only an approximation of the true derivative of the cost function when using Glen's flow law (Eq.1) with n > 1 in Eq. (2).

Regularisation
To avoid non-physical negative values of the basal friction coefficient, β is expressed as The optimisation is now done with respect to α and the Gâteaux derivative of J o with respect to α is written as To avoid over fitting of the data and improve the conditioning of the problem, a smoothness constraint is added to the cost function under the form of a Tikhonov regularisation term penalising the spatial first derivatives of α as The Gâteaux derivative of J reg with respect to α for a perturbation α is obtained by The total cost function is now written as where λ is a positive ad-hoc parameter.The minimum of this cost function is no longer the best fit to observations, but a compromise (through the tuning of λ) between fit to observations and smoothness in α.

Minimisation
Surface ice flow speeds vary over several order of magnitudes between the interior of the ice sheet (slow flow) and the outlet glaciers (fast flow).Therefore, Schäfer et al. (2012) have shown that good convergence properties are obtained with the Robin inverse method by using a spatially varying step size rather than the fixed-step used in the original gradient descent algorithm of Arthern and Gudmundsson (2010).
We apply here another strategy where the Gâteaux derivatives given by a continuous scalar product represented by the integral on b , Eqs. ( 17) and ( 19), are transformed in a discrete euclidean product when discretized on the finite element mesh.The area surrounding each finite element node on b is then included in the gradients used for the minimisation.This leads to good convergence properties with our unstructured mesh as large elements correspond to low velocity areas, and vice versa.The minimisation of the cost function J tot with respect to α is done using the limited memory quasi-Newton routine M1QN3 (Gilbert and Lemaréchal, 1989) implemented in Elmer in reverse communication mode.This method uses an approximation of the second derivatives of the cost function and is then more efficient than a fixed-step gradient descent.

Results
The observed velocities, shown in Fig. 1a, are a compilation of data sets obtained from RADARSAT data at different dates during the first decade of the 2000s (Joughin et al., 2010).We choose this compilation as it gives the best coverage.For glaciers that have been accelerating, it therefore provides a kind of average value for this period.But for these glaciers, the surface topography has also diverged from the surface topography used here.A proper comparison of the model results with observations of the 2000s would therefore require coherent data sets for both the topography and the velocities, which are currently not available.For the inversion, gaps in data have been filled by the balance velocities available in the SeaRise data set.As these gaps are mainly located in areas of low speed in the central parts of the ice sheet, except in the south and south-east, no special care has been taken to insure continuity between balance and observed velocities.
To compare the performance of the two inverse methods and assess the dependence of the results to the initialisation of the basal friction coefficient β, the initial field is given by for the Robin inverse method where U bal is the balance velocity, and by for the control inverse method.
The optimal regularisation parameter λ in Eq. ( 20) is chosen using the L-curve method (Hansen, 2001).The L-curve is a plot of the optimised variable smoothness, i.e. the term J reg in Eq. ( 20), as a function of the mismatch between the model and the observations, i.e. the term J o in Eq. ( 20).The L-curve obtained with both methods is given in Fig. 3.The term J reg is identical for the two methods, whereas J o is not (Eqs.12 and 14).This leads to different values of the regularisation parameter λ in Eq. ( 20), which allows tuning the weight of the regularisation with respect to J o .
For the Robin inverse method, when increasing λ from 0 to 10 9 , the roughness of the basal friction coefficient field, represented by J reg , decreases by several orders of magnitude with a small decrease of the mismatch between the model and the observation, represented by J o .This decrease of J o may be due to the fact that the gradient used in the model, Eq. ( 13), is not the true gradient of J o due to the nonlinearity of the ice rheology.The regularisation, for which the gradient is exact, therefore improves the convergence properties of the model.For higher values of λ, J reg still decreases but with a concomitant increase of J o as the basal friction coefficient field becomes too smooth.The L-curves obtained with the two methods are very similar.The minimum J o is obtained for the same order of magnitude of J reg .The optimal value for λ is chosen as the minimum value of J o , i.e. λ Robin = 10 8 for the Robin inverse method, and λ CI = 10 11 for the control inverse method.
The surface velocity field obtained after optimisation of the basal friction coefficient field with the Robin inverse method is shown in Fig. 2. Our implementation reproduces very well the observed large-scale flow features (Fig. 1a) with low velocities in the interior and areas of rapid ice flow, restricted to the observed outlet glaciers, near the margins.The largest outlet glaciers (Jakobshavn Isbrae, Kangerlugssuaq, Helheim, . . . ) and their catchments are well captured by the anisotropic mesh and the modelled velocity pattern is in good agreement with the observations.Smaller outlet glaciers down to a few kilometres in width are also individually distinguishable.
The absolute and relative errors on the surface velocities at the end of the optimisation are shown in Figs. 4 and 5, respectively.Similar to the velocity magnitude, the absolute error varies by several orders of magnitude between the interior and the margins.The relative error is only few percents in most of the interior where ice is flowing faster than few meters per year.This error is usually higher very locally near the margins.The highest relative errors are located in the North in Petermann Glacier and in the Northeast Greenland Ice Stream where long floating tongues are present but not explicitly taken into account in this application of the model.The remaining differences between modelled and observed velocities can come from four main reasons: i. non convergence of the minimisation: it has been shown on twin experiments where the minimum is known (Arthern and Gudmundsson, 2010) that the gradients of the cost function derived analytically for a linear rheology, Eqs. ( 13) and ( 15), work well in practice with a non-linear rheology.But there is no guarantee that the actual minimum will be found (Goldberg and Sergienko, 2011), especially in real applications where the curvature of the cost function is very low.
ii. Remaining uncertainties: adjusting the sliding coefficient can compensate only partly for errors associated with the uncertainties on the other model initial conditions.iii.Too coarse spatial resolution of the model where the minimum mesh resolution is lower than those of the velocity data, so that the model will not be able to capture all details, especially for the smallest outlets; iv.Too coarse spatial resolution of the data as, for example, the ice thickness is not sufficiently known in most outlet glaciers.
It is difficult to test the first hypothesis as the minimum is unknown but both the cost function and the norm of the gradient decrease during the minimisation and both inverse methods lead to very similar results (not shown for the control inverse methods), so that we are confident to be close to the actual minimum.
Errors shown in Figs. 4 and 5 account both for the error on the direction and on the magnitude of the modelled velocities compared to the observations.The direction of the flow is mainly governed by the ice-sheet topography and adjusting the sliding coefficient has little effect on the flow direction.The cost function used with the control inverse method only accounts for the difference between the velocity norms and not for the direction.Both inverse methods lead to similar errors so that we assume that most of the absolute error is representative of the error on the velocity norm.For the outlets, a common feature is that model velocities are lower than the observations along the central flow lines but higher along the shear margins (not shown here).This can be explained by insufficient resolution of the model and/or of the data, or remaining uncertainties on the ice viscosity for example.

Relaxation
When running prognostic simulations from this initialisation, the free-surface elevation shows non-physical rates of change especially at the margins (Fig. 6).The free surface of the ice sheet is therefore allowed to relax during a 50-yr timedependent run, forced by a constant present-day surface mass balance field (Ettema et al., 2009).
One contribution to these initial surface changes comes from local ice flux divergence anomalies (Seroussi et al., 2011) due to uncertainties on the model initial conditions, including initial topography and model parameters.These anomalies disappear very quickly within a few years.The other contribution is a response to the fact that many outlet www.the-cryosphere.net/6/1561/2012/The Cryosphere, 6, 1561-1576, 2012 glaciers are not able to evacuate the ice flowing through them and they need to thicken up.One explanation is that the ice thicknesses of many outlet glaciers is poorly known and they may be thicker in reality than the current geometry data allows for.The other explanation is that the lateral side of the mesh does not correspond exactly to the current position of the fronts of marine terminating glaciers for which we have no compilation for the whole GrIS.The ice then needs to be transported towards the edge of the domain to be evacuated as discharge through the lateral boundary.This effect can be shown on the results for the ice discharge globally (Fig. 7) and locally for individual outlet glaciers (Table 2).Ice discharge is very low and underestimated for most of the smallest outlet glaciers at the beginning of the relaxation but increases with time as ice reaches the edges of the domain.
The time needed to offset these uncertainties and reach an equilibrium is difficult to estimate in advance and will be different for each drainage basin.The relaxation duration is arbitrarily fixed to 50 years.At the end of the relaxation, the surface velocity structure, given in Fig. 1b, remains similar to the observations, with rapid ice flow areas still concentrated in known outlet glaciers.The various terms of the mass balance equation during the relaxation are given in Fig. 7. SMB is nearly constant and equal to 430 Gt a −1 during the relaxation showing that the glaciated area has not changed dramatically.Most of the ice discharge increase happen during the first twenty years to reach 200 Gt a −1 , after this increase slows down.During this relaxation the ice volume increases by less than 0.5 %.Running the model longer should lead to a steady state where the ice discharge should offset the sur-Table 2. Model discharge for individual outlets and the whole ice sheet in Gt a −1 : after 1 yr (RI 1a) and after 50 yr (RI 50a) of surface relaxation with the Robin inverse method, after 50 yr of surface relaxation with the control inverse method (CI 50a), and observations from 1996 (or 2000 when not available, and converted from km 3 a −1 to Gt a −1 using a uniform density of 910 kg m −3 for ice) (Rignot and Kanagaratnam, 2006)  face mass balance as shown in Fig. 8 when running the same constant climatic scenario for another century.
At the end of the relaxation, the rate of increase of the ice discharge is around 1 Gt a −2 and the free surface is nearly at equilibrium except in a few areas near the margins.As shown by Fig. 9a when running the same scenario for another 10 years, this remaining elevation change is mainly positive as these areas are still unable to evacuate all the ice flowing to them but the rate of change is within a physically acceptable range (Fig. 6).The change in surface elevation between the beginning and the end of the relaxation exceeds several hundred meters in some places.This difference is large but is of the same order of magnitude as the uncertainty on the ice thickness in some areas of the GrIS margins.This remains the main limitation of the model where only the upper surface is adjusted to compensate the remaining uncertainties in the model initial conditions.
During this relaxation period with a constant forcing, each drainage basin has been driven towards a steady state and the model results should preferentially be compared with observations before the recent glacier accelerations.Computed discharge values from the main outlet glaciers are compared with observations given by Rignot and Kanagaratnam (2006) for 1996, or 2000 when not available.A detailed comparison shows that the modelled ice discharge is generally underestimated.For the three largest outlet glaciers, where the topography is the best known, the agreement is within 1 % for Helheim Glacier and 15 % for Jakobshavn Isbrae, but D is largely underestimated by more than 50 % at Kandgerdlugssuaq where the drainage basin still shows high rates of surface-elevation change (Fig. 6).Discharge at Petermann Glacier and the Northeast Greenland Ice Stream is also largely underestimated by up to 80 % and these glaciers are still thickening at high rates.This can be explained by the fact that their floating tongues are not treated explicitly in the model resulting in the neglect of the mass loss from melting below the shelves at the interface with the ocean.The total computed discharge for the relaxed solution is around 300 Gt a −1 .Up to now, the magnitude of the computed ice discharge in a continental ice-sheet model has only been addressed by Price et al. (2011) by tuning the boundary condition at the ice front to reproduce the observations only on three major outlets.The evolution of the total volume and ice discharge obtained with the basal friction coefficient field optimised with the two inverse methods are very similar.This is true for each individual outlet as shown in Table 2.The two inverse methods perform similarly and neither can be favoured in view of these results or in terms of computation performance.

Sensitivity experiments
We use the relaxed solution of the Robin inverse method as the starting point to investigate the GrIS mass balance over one century.

Set up
The climate forcing in the model impacts ice dynamics only through the net accumulation rate at each surface mesh node as the ice temperature field is kept constant.Two SMB scenarios are used, the first (C1), corresponds to keeping present conditions (Ettema et al., 2009) constant with time; the second (C2) represents an ensemble of 18 climate models forced under the IPCC A1B emission scenario.These forcings are taken from the freely available SeaRISE data compilations (http://tinyurl.com/srise-umt).
Here, we especially focus on the GrIS dynamical response to an increase in basal lubrication.Perturbations are introduced by applying three homogeneous changes in the basal friction coefficient β from the initial field β i (x, y) inferred by the inverse method (Sec.3.1).
i.No perturbation with β unchanged form the initial field (BF1): ii.A constant perturbation, introduced as a step function, with β divided by two and then kept constant (BF2): iii.A continuously enhanced perturbation with β reduced by one order of magnitude over one century (BF3): where t is the simulation time (0 ≤ t ≤ 100) and β = 10 α .As has been shown for a surging glacier (Jay-Allemand et al., 2011), β can vary by several orders of magnitude with only small changes in the water pressure so that the magnitude of the imposed changed would not be unrealistic if applied in few locations and for short time periods.However, it seems highly unlikely that the GrIS will surge all at the same time, so that this last sensitivity experiment is used as an high end scenario of the possible future for the GrIS mass loss.
Further destabilisations introduced by changes in the seaward boundary conditions or weakening of the lateral margins are excluded from this study.Experiments are hereafter referred to by the climate forcing (C1 or C2) and the basal friction scenario (BF1, BF2 or BF3).

Results
We evaluate the results of the perturbations by considering ice-flow velocities (Fig. 1), the ice-sheet total mass balance (Fig. 8), discharge values from the main outlet glaciers (Table 3) and free-surface-elevation rate-of-change (Fig. 9).For the ice-sheet total mass balance (Fig. 8), we present the total ice-sheet volume and volume change from the starting point, converted to sea-level equivalents (Fig. 8a).The annual mass balance (MB, Fig. 8b) is obtained as the time derivative of the ice-sheet volume, and the respective contributions of D and SMB are given in Fig. 8c.The difference between MB and SMB-D is the ice flux lost through the bottom boundary (see Sect. 2.1, discussion on boundary condition Eq. 7).In all the applications it corresponds approximately to 10 % of D.
With constant conditions (C1 BF1), the model tends to reach a steady state where the modelled discharge balances the current SMB (430 Gt a −1 , Fig. 8c).Neither the glaciated area, nor the surface velocity pattern, change dramatically during this experiment (Fig. 1c) and the ice discharge shows a small increase in the main outlets (Table 3).
The climatic perturbation used here (SMB scenario C2) shows a reduction of SMB of approximately only 100 Gt a −1 after one century (Fig. 8c), which is a lower bound of the forecast given by current climate models (Fettweis et al., 2008).Changes in the glaciated area between the three perturbation experiments detailed below lead to only small differences in the total SMB after one century (Fig. 8c).These differences are one order of magnitude lower than those between the various climate models (Fettweis et al., 2008) but do not account for the possible feedback between the icesheet topography and climate as the prescribed accumulation is only a function of position and time.Investigating this feedback requires proper coupling of high resolution icesheet and climate models.
Without dynamic perturbation (C2 BF1), the computed ice discharge is of the same order of magnitude as the SMB and decreases at an equivalent rate (Fig. 8c).The resulting annual mass balance is nearly centred around zero and lacks a particular trend over the century (Fig. 8b), as a consequence, the total ice volume is nearly constant (Fig. 8a).For this experiment, the velocity pattern (Fig. 1d) remains similar to the present one.The increase of ablation for this climate scenario is higher in West Greenland.This leads to thinning of the marginal ice in this area, resulting in a retreat of land terminating glaciers and in a decrease in ice discharge of marine terminating glaciers (Table 3).With a constant dynamic perturbation (C2 BF2), halving β before the simulation results in an almost immediate doubling of the ice discharge bringing it close to 500 Gt a −1 (Fig. 8c).When looking at the total discharge only, this perturbation allows to compensate the model initial uncertainties as this value is more in agreement with current estimates based on observations (Rignot et al., 2011).However, this is done at the expense of the velocity pattern as we see higher velocities on the interior of each drainage basin (Fig. 1e) compared to the current observations (Fig. 1a).A progression of high velocitie areas upstream of each drainage basin may be expected as ablation will reach higher elevations in the future, increasing water pressure and basal lubrications in these areas.However, the expected effect in low elevation areas is not clear as the increase of water availability should contribute in the formation of efficient drainage systems resulting in a decrease of the water pressure and basal lubrication (Schoof, 2010).With this experiment, we investigate how an initial ice-sheet imbalance of −200 Gt a −1 , which is close to the current observations, can evolve within a century.The computed total discharge decreases throughout the 100yr simulation at a rate equivalent to the rate of decrease of SMB except during the first decade where it decreases faster probably as a reaction to the initial perturbation (Fig. 8c).As a result, after the first decade the annual mass balance shows no particular trend (Fig. 8b) and the ice sheet loses mass at a nearly constant rate throughout the century (Fig. 8a).With this experiment, the retreat of land terminating glaciers in the west coast induced by the increase of ablation is compensated by their acceleration, allowing to drain more ice to the areas subjected to high ablation.Again, the reduction of the discharge throughout the simulation is higher in the west coast as more ice is melted before reaching the margins of the domain.
With an increasing dynamic perturbation (C2 BF3), the discharge increases continuously from 300 Gt a −1 to 1400 Gt a −1 after 100 yr (Fig. 8c), and the ice sheet is losing mass throughout the 100-yr simulation (Fig. 8a, b).With this result we show that the acceleration of the ice discharge observed during the last decade (Rignot et al., 2011) can be maintained over a whole century.However, the velocity pattern at the end of the century (Fig. 1f), shows large areas of high velocities (>100 m a −1 ) far inland of each glacier.Such a scenario seems unlikely as it would required a sustained surge of the whole ice sheet in the same time, which is excluded if efficient drainage systems develop as a result of increased water supply (Schoof, 2010).As a consequence, we estimate that, in the absence of further destabilisation of the GrIS ocean margins, this scenario represents an upper bound for the ice-sheet contribution to SLR by 2100.
An important output of these experiments is that the ice discharge shows an extremely rapid adjustment to perturbations as shown at the beginning of the relaxation (Fig. 7b) and in response to the dynamical perturbations (for both experiments BF2 and BF3, Fig. 8c).However, the spatial and temporal variability of the forcings affecting the basal and seaward boundary conditions of the ice sheet remain poorly constrained, limiting the forecast ability of the models.
Together with surface velocities, rates of surface elevation change are now widely measured and available (Pritchard et al., 2009).All the glaciers that have been accelerating recently also show a dynamical thinning.Even if these observations have been used by flow models as a post-validation to try to discriminate between the destabilising processes (Nick et al., 2009;Price et al., 2011), they have not been used in a proper inverse method so far.Because the homogeneous dynamical perturbation applied here is probably too crude and does not account for the observed retreat of marine terminating glaciers, our results are not directly compared with observations.However, we provide a qualitative discussion of the surface elevation changes computed after 10 yr of simulations (Fig. 9).With experiment C1 BF1, as discussed previously, some drainage basins are still not at equilibrium and their outlets are still thickening due to a SMB higher than the computed discharge.In the interior, the surface elevation change is nearly zero.With a climate perturbation only, i.e. experiment C2 BF1, the margins in the west and south-east are thinning and the thickening of Petermann Glacier and the Northeast Greenland Ice Stream is less pronounced.These differences with experiment C1 BF1 mostly come from the change in the SMB term.
Experiments with a dynamical perturbation, i.e C2 BF2 and C2 BF3, show an additional dynamical thinning associated with the acceleration of the ice, upstream of each drainage basin.Downstream, near the margins, this gives two different behaviours.(i) If the decrease of the basal friction coefficient also produces an acceleration and an increase of the discharge sufficient to offset the mass excess coming from upstream, the whole ice stream shows a dynamical thinning.This is the case in the south-east and for Jakobshavn Isbrae and Heilhem Glacier.(ii) If the acceleration and the increase of ice discharge are not sufficient, this results in a dynamical thickening of the margins.This is the case in the north-west and for Kangerdlussuaq Glacier.For land terminating glaciers south of Jakobsahvn Isbrae, the thinning induced by the increase of ablation is partly offset by this dynamical thickening.
In some aspects, these results therefore agree with observations and show an interesting regional variability of the response to the dynamical perturbation.The link between surface runoff, basal hydrology and basal sliding will have to be investigated more in-depth to make quantitative comparisons.

Conclusions
We have shown that our implementation of a model with the correct treatment of longitudinal stresses, sufficient resolution to resolve medium to small outlet glaciers, and a careful initialisation of ice-flow parameters allows us to satisfy the essential prerequisite of simulating the observed velocity structure of the ice sheet.Neither of the two recent inverse methods used to infer the basal friction coefficient field from the knowledge of the ice-sheet topography and ice-surface velocities can be favoured in view of our results.
Due to the remaining uncertainties in the model initial conditions, the ice-sheet surface has been allowed to diverge from the observations during a relaxation period of 50 years.This is the main limitation of the model as this relaxation drives the model towards a steady state, limiting our ability to capture the currently observed transient response of the GrIS.However, we have shown that, after the relaxation, the ice discharge obtained globally and for individual outlet glaciers is still in good agreement with observations of the late 1990s.The main remaining uncertainties in the model are the ice viscosity field, the bedrock elevation and the position of the ice margins.Future work will benefit from new observations and from the developments of inverse methods to constrain these uncertainties (Arthern and Gudmundsson, 2010;Pralong and Gudmundsson, 2011;Morlighem et al., 2010;Petra et al., 2012).However, these methods remain restricted to diagnostic simulations and proper methods, such as automatic differentiation of the true discrete adjoint model (Heimbach and Bugnion, 2009), still have to be applied to higher-order ice flow modelling to capture the observed transient response of the ice sheet.
More specific projections will arise as the ice sheet is driven by more complete and precise climate scenarios, and with greater understanding of processes at the front of outlet glaciers.Our new-generation continental-scale ice-sheet model is well suited to incorporate such information as it becomes available.However, our current model experiments are significant and general conclusions can already be drawn.
Results confirms that the overall mass balance of the GrIS is sensitive, not only to changing SMB, but in large parts, to ice discharge.Our model shows a rapid adjustment of the ice discharge in response to dynamical perturbations.This is supported by other models that include processes at the ice front (Nick et al., 2009(Nick et al., , 2012) ) and by recent observations (Howat et al., 2007).Moving ice margins and a physically based treatment of calving are essential future developments of the model, and are required to move from sensitivity experiments to realistic projections of the ice discharge evolution.
Results show that, unless the perturbation is continuously enhanced, an increase of the surface ablation reduces the discharge, stabilising the ice sheet.However, this effect could be negated if the thinning of outlet glaciers results in further destabilisation and retreat of the ocean margins making it possible to sustain high discharge rates.Nevertheless, this demonstrates that, even on a century timescale, discharge and surface mass balance anomalies cannot be treated independently in estimating the GrIS contribution to SLR.This effect could be even higher due to the possible feedback between surface elevation and surface mass balance.Proper coupling with local surface mass balance models will then be required to improve the model predictive ability.
We use the modelling conducted here and the currently observed ice-sheet mass balance to estimate possible bounds for the future GrIS contribution to SLR.These bounds do not take into account possible further movements of the GrIS ocean margins as proper treatment of the processes affecting their stability in a whole ice-sheet model is still lacking.
In our experiments with constant basal conditions, the rate of decrease of D and SMB are similar, resulting in a nearly constant annual mass balance.In these conditions, a probable lower bound for the ice-sheet contribution to SLR is given by assuming that the currently observed ice-sheet imbalance will be preserved over the whole century.Taking an averaged value of −300 Gt a −1 for the ice-sheet mass The Cryosphere, 6, 1561-1576, 2012 www.the-cryosphere.net/6/1561/2012/balance in 2010 (Rignot et al., 2011) leads to a total contribution to SLR of 75 mm by 2100.This value is comparable to the 46 mm (40 mm from SMB and 6 mm from dynamics) given by Price et al. (2011) as an estimate of the committed SLR due to the observed last decade ice-sheet dynamics.Our value is closer to their estimated upper bound (85 mm including 45 mm from dynamics) which is estimated by assuming a self-similar response of the ice sheet to a 10 yr-recurring forcing in the future decades (Price et al., 2011) .
We have shown that, conversely, if the forcing is continuously increased, there is sufficient ice available to sustain the current rate of increase in discharge over an entire century.However, in our model, this is achieved by a very strong perturbation of the basal lubrication field.Under these conditions, a probable upper bound for the ice-sheet contribution to SLR is given by assuming that this rate of increase will be preserved in the future.Taking a discharge anomaly of 100 Gt a −1 increasing at rate of 9 Gt a −2 for year 2010 (Rignot et al., 2011), leads to a contribution to SLR of 140 mm by 2100 from dynamics only.This value is in the lower half of the values obtained from kinematic considerations assuming low (93 mm) and high (467 mm) scenarios (Pfeffer et al., 2008).

Fig. 2 .
Fig. 2. Unstructured finite element mesh and model surface velocities after optimisation of the basal friction coefficient with the Robin inverse method.Colored boxes show close-up views for various outlet glaciers of interest.

Fig. 4 .
Fig. 4. Absolute error on surface velocities |u mod − u obs | in m a −1 at the end of the optimisation using the Robin inverse method.

Fig. 5 .
Fig. 5. Relative error on surface velocities |u mod − u obs |/|u obs | in % at the end of the optimisation using the Robin inverse method.Areas where |u obs | < 2.5 m a −1 have been removed from display.

Fig. 6 .
Fig. 6.Free surface elevation absolute rate of change: (a) after 1 yr of relaxation; (b) after 50 yr of relaxation.

Fig. 7 .
Fig. 7. GrIS mass balance during relaxation: (top) evolution of the total ice volume in meters of sea level equivalent and (bottom) evolution of the discharge (solid lines), SMB (dashed lines) and flux through the bedrock (dotted lines) in Gt a −1 .

Fig. 8 .
Fig. 8. GrIS future mass balance as a function of time: (a) total ice volume (left axis, meters of sea level equivalent) or volume change from initial time (right axis, milimeters of sea level equivalent), (b) annual mass balance (the brown line correspond to 0) and (c) discharge (D) (solid lines) and SMB (dashed lines) in Gt a −1 .

Table 1 .
List of parameter values used in this study. (Obs.).

Table 3 .
Model discharge for individual outlets and the whole ice sheet in Gt a −1 after one century for the various experiments.