Simulations of the Greenland ice sheet 100 years into the future with the full Stokes model Elmer/Ice

Abstract

It is likely that climate change will have a significant impact on the mass balance of the Greenland ice sheet, contributing to future sea-level rise. Here we present the implementation of the full Stokes model Elmer/Ice for the Greenland ice sheet, which includes a mesh refinement technique in order to resolve fast-flowing ice streams and outlet glaciers. We discuss simulations 100 years into the future, forced by scenarios defined by the SeaRISE (Sea-level Response to Ice Sheet Evolution) community effort. For comparison, the same experiments are also run with the shallow-ice model SICOPOLIS (SImulation COde for POLythermal Ice Sheets). We find that Elmer/Ice is ~43% more sensitive (exhibits a larger loss of ice-sheet volume relative to the control run) than SICOPOLIS for the ice-dynamic scenario (doubled basal sliding), but ~61 % less sensitive for the direct global warming scenario (based on the A1 B moderate-emission scenario for greenhouse gases). The scenario with combined A1B global warming and doubled basal sliding forcing produces a Greenland contribution to sea-level rise of ~15cm for Elmer/Ice and ~12cm for SICOPOLIS over the next 100 years.

1. Introduction

The Greenland ice sheet is the second largest land ice mass on the present-day Earth, and its volume amounts to ~7.3 m s.l.e. (metres sea level equivalent). The current mass balance of the ice sheet is most likely negative with an accelerating trend, though the uncertainty is significant (Lemke and others, 2007; Rignot and others, 2011). Surface melting increases strongly with rising surface temperatures, making the ice sheet very susceptible to future global warming. In addition, recent observations (Zwally and others, 2002; Rignot and Kanagaratnam, 2006; Howat and others, 2007; Joughin and others, 2008) have led to strong concerns that ice-dynamical processes (basal sliding accelerated by surface meltwater, speed-up of ice streams and outlet glaciers) may boost the decay and thus lead to an additional contribution to sea-level rise. Therefore, it is clearly necessary to comprehensively model the dynamics of the Greenland ice sheet, including ice streams and outlet glaciers.

Many models have been developed to simulate the dynamics and evolution of ice sheets and glaciers. The shallow-ice approximation (Hutter, 1983; Morland, 1984) has been widely used for ice-sheet models (e.g. Huybrechts, 1990; Calov and Hutter, 1996; Ritz and others, 2001; Saito and Abe-Ouchi, 2004; Rutt and others, 2009). This approximation neglects the normal deviatoric stress and horizontal shear components, thus implying a significant simplification that works well for large-scale ice-sheet dynamics but is inappropriate in the vicinity of ice divides and margins, fast-flowing regions like ice streams, and small steeply inclined glaciers in general (e.g. Greve and Blatter, 2009). This gave rise to the formulation of higher-order models (Blatter, 1995; Baral and others, 2001; Hindmarsh, 2004) in which longitudinal stresses are taken into account to various extents. Many of these models have been applied to two-dimensional (2-D) domains (Dahl-Jensen, 1989; Blatter, 1995; Colinge and Blatter, 1998; Pattyn, 2000), and Dahl-Jensen (1989) demonstrated the importance of longitudinal deviatoric stresses for plane flow along a flowline. Pattyn (1996, 2000) and Pattyn and Decleir (1998) applied a 2-D higher-order model with thermomechanical coupling to Shirase drainage basin, Dronning Maud Land, Antarctica. Pattyn (2003) developed a three-dimensional (3-D) higher-order thermomechanical sheet model, carried out the European Ice-Sheet Modelling Initiative (EISMINT) I and II benchmark experiments (Huybrechts and others, 1996; Payne and others, 2000) and provided a comparison with the Saito-Blatter model that also includes higher-order dynamics (Saito and others, 2003). More recent developments are the models of Pollard and DeConto (2007, 2009) and Bueler and Brown (2009) that employ heuristic combinations of the shallow-ice and shallow-shelf approximations (Morland, 1987; MacAyeal, 1989), as well as the application of a first- order model to the Greenland ice sheet by Price and others (2011). However, due to several shortcomings inherent in those models, none of their results contributed to the Fourth Assessment Report (AR4) of the Intergovernmental Panel on Climate Change (Solomon and others, 2007), which represents a great opportunity for the development and application of full Stokes models.

Models that solve the full Stokes equations (in which all stress components are accounted for) in two or three dimensions have been proposed and applied mainly to glacier systems or parts of an ice sheet (e.g. Gudmundsson, 1999; Sugiyama and others, 2003; Martfn and others, 2004; Price and others, 2007; Jouvet and others, 2009). Comparisons between various full Stokes and higher-order models were carried out in the Higher-Order Model (HOM) intercomparison topic of the Ice-Sheet Model Intercomparison Project (ISMIP) (Pattyn and others, 2008). It was found that all participating models produced results that are in close agreement. However, the full Stokes models were most consistent with each other, whereas the spread among the various higher-order models was larger, thus clearly motivating the use of full Stokes models.

Apart from the recent studies of Ren and Leslie (2011) and Ren and others (2011a, b), full Stokes models have not yet been applied to an entire ice sheet because of the enormous computational demand. Here the full Stokes thermomechanically coupled model Elmer/Ice (e.g. Zwinger and others, 2007; Gagliardini and Zwinger, 2008; Durand and others, 2009; Zwinger and Moore, 2009; Seddik and others, 2011) is applied to the Greenland ice sheet. Elmer/Ice employs the finite-element method to solve the full Stokes equations, the temperature evolution equation and the evolution equation of the free surface. The general framework of this modelling effort is a contribution to the Sea-level Response to Ice Sheet Evolution (SeaRISE) assessment project, a community-organized effort to estimate the likely range of ice-sheet contributions to sea-level rise over the nextfew hundred years (http://tinyurl.com/srise-lanl, http://tinyurl.com/srise-umt). We therefore carry out the four SeaRISE experiments considered by Greve and others (2011), who defined climatic and dynamic future scenarios. Results are also compared with the shallow-ice approximation model SICOpOlIS (SImulation COde for POLythermal Ice Sheets (e.g. Greve, 1997, 2000; Greve and others, 2011)) in order to assess the differences in the response of the two models.

2. ELMER/ICE: Thermomechanically Coupled Full Stokes Flow Model

2.1. Dynamic/thermodynamic model equations

2.1.1. Field equations

Since ice is an (almost) incompressible material, conservation of mass requires that the velocity field (vector v) is solenoidal,

(1)

Further, the acceleration (inertia force) is negligible, so the equation of motion is given by the incompressible Stokes equation,

(2)

(e.g. Greve and Blatter, 2009), where p is the pressure, n the viscosity, ρ the ice density and g =−ge_z the gravitational acceleration vector pointing downward. The viscosity is described by Glen’s flow law,

(3)

where is the effective strain rate, D = sym the strain-rate tensor (symmetric part of the velocity gradient L = grad v), n the power-law exponent, T’ = T− T_m the temperature relative to pressure melting (T is the absolute temperature, T_m = T₀− βp is the pressure melting point, T_o is the melting point at low pressure and j3 is the Clausius-Clapeyron constant), A(T’) the rate factor and E the flow-enhancement factor. The rate factor is expressed by the Arrhenius law

(4)

where A ₀ is the pre-exponential constant, Q the activation energy and R the universal gas constant. All parameters are given in Table 1.

Table 1. Standard physical parameters used for the simulations with both Elmer/Ice and SICOPOLIS

The temperature equation follows from the general balance equation of internal energy and reads

(5)

(e.g. Greve and Blatter, 2009), where k and c are the heat conductivity and specific heat of ice, respectively (Table 1).

The free surface equation follows from the kinematic boundary condition formulation and reads

(6)

where h(x, y, t) is the free surface and a_s(x, y, t) is the accumulation-ablation function or surface mass balance. The ice base, b(x, y), is assumed to be rigid (isostatic compensation neglected) and thus at all times equal to the prescribed initial condition.

2.1.2. Boundary conditions

We extract the boundary conditions required to close the system of equations posed in Section 2.1.1 mainly from the SeaRISE specifications (see also Greve and others, 2011). The ice surface is assumed to be stress-free (atmospheric pressure and wind stress neglected). The surface air temperature is parameterized as a function of surface elevation, h, latitude ϕ, longitude, λ, and time, t, following Fausto and others (2009):

(7)

where T _ma and T _mj are the mean annual and mean July (summer) surface temperatures, respectively, the temperature constants are d _ma = 41.83°C and d _mj = 14.70°C, the mean slope lapse rates are γ_ma = −6.309°Ckm^-1 and γ_mj = −5.426°Ckm^-1, the latitude coefficients are c _ma = −0.7189°C(°N)^-1 and c _mj = −0.1585°C(°N)^-1, and the longitude coefficients are K _ma = 0.0672°C(°Wr¹ and K _mj = 0.0518°C(°W)^-1.

The purely time-dependent anomaly term, ΔT(t), describes the deviation from present-day conditions. For the past, it is based on the oxygen isotope record (δ¹⁸O) from the Greenland Icecore Project (GRIP) ice core (Dansgaard and others, 1993; Johnsen and others, 1997), which was converted to a record of temperature variation from 125 ka BP to the present (here, the notation ka BP means thousand calendar years before present). For the future, Eqn (7) is only used for the experiments with constant, present-day climate forcing (thus ΔT(t) = 0), whereas the experiments with AR4 climate forcing are driven directly by an ensemble average of simulated surface temperatures (Section 3.2).

For the present-day mean annual precipitation rate, P _ma,present(λ, ϕ), recent data of Ettema and others (2009) are used. Past precipitation rates are not required in this study because of the fixed-topography spin-up approach (Section 3.1). For the future runs with constant, present-day climate forcing, P _ma,present(λ, ϕ) is used unchanged, while the AR4 climate experiments are driven directly by simulated precipitation rates, analogous to the surface temperature.

Surface melting is parameterized by Reeh’s (1991) positive degree-day (PDD) method, supplemented by the semianalytical solution for the PDD integral by Calov and Greve (2005). The PDD factors are β_ice = 8 mm (ice) d^-1 °C^-1 for ice melt and β_snow = 3 mm (ice) d^-1 °C^-1 for snowmelt (Huybrechts and de Wolde, 1999). Furthermore, the standard deviation of short-term, statistical air-temperature fluctuations is σ = 5°C (Huybrechts and de Wolde, 1999), and the saturation factor for the formation of superimposed ice is chosen as P _max = 0.6 (Reeh, 1991). Conversion from the present-day mean annual precipitation (Ettema and others, 2009) to the snowfall rate (solid precipitation) is done on a monthly basis using the empirical relation of Marsiat (1994). Mean monthly rainfall (liquid precipitation) is obtained as the difference between precipitation and snowfall.

At the base, described by the function z = b(x, y), a Weertman-type sliding law with sub-melt sliding is used (Greve, 2005),

(8)

where τ_b is the basal drag (shear stress), N _b the basal normal stress, T’_b the basal temperatures relative to pressure melting, is the sliding coefficient, p = 3, q = 2 are the sliding exponents and γ = 1°C is the sub-melt-sliding parameter (Hindmarsh and Le Meur, 2001).

If the base is cold (temperature below the pressure-melting point), the energy jump condition yields a Neumann-type boundary condition for the basal temperature,

(9)

(e.g. Greve and Blatter, 2009), where the geothermal flux distribution over the bedrock is given by Shapiro and Ritzwoller (2004) and t is the Cauchy stress tensor. If the base is temperate (temperature at the pressure-melting point), the energy jump condition determines the basal melting rate,

(10)

where L is the latent heat of ice.

At the lateral boundaries of the domain (vertical faces at the ice margin), the stress-free condition is applied, and the horizontal temperature gradient is assumed to vanish (zero- flux condition).

2.2. Finite-element implementation

The model equations detailed in Section 2.1 are solved numerically with the Elmer/Ice model. It is based on the open-source multi-physics package Elmer developed at the CSC− IT Center for Science in Espoo, Finland (http://www.csc.fi/elmer/), and uses the finite-element method.

The model domain covers the entire area of the present- day Greenland ice sheet. The domain is projected to a polar stereographic map with standard parallel 71° N and central meridian 39 ^° W. The present geometry (surface and basal topographies) is derived from U.C. Herzfeld and others (unpublished information) where the basal topography was created so that the troughs atJakobshavn Isbr? and Helheim, Kangerdlugssuaq and Petermann glaciers are preserved. A mesh of the computational domain is created using an initial footprint that contains elements of 5 km horizontal resolution. To limit the number of elements on the footprint while maximizing the spatial resolution in regions where physics demands higher accuracy, an anisotropic mesh adaptation scheme is employed. Its metric is based on the Hessian matrix of the observed surface velocities (distributed by SeaRISE based on work by Joughin and others (2010) with gaps filled by balance velocities of Bamber and others (2001b)) in order to equi-distribute the a priori error estimate using an edge-based anisotropic mesh optimization. The metric tensor is computed following Frey and Alauzet (2005), and the adaptation is carried out with the automatic tool YAMS (Frey, 2001; Morlighem and others, 2010). The resulting mesh in the central part of the Greenland ice sheet, including the refinements at Jakobshavn Isbr* and Kangerdlugssuaq Glacier, is depicted in Figure 1. The final footprint is vertically extruded to form a 3-D mesh of 320 880 elements with 17 equidistant, terrainfollowing layers.

Fig. 1. Observed surface velocities of the central part of the Greenland ice sheet (distributed by SeaRISE based on work by Joughin and others (2010) with gaps filled by balance velocities of Bamber and others (2001b)) and anisotropic mesh with the clearly visible refinements at Jakobshavn Isbræ (JIS) and Kangerdlugssuaq Glacier (KL).

The nonlinearity of the model equations is dealt with by a Picard iteration scheme, and stabilization methods (Franca and Frey, 1992; Franca and others, 1992) are applied to the finite-element discretization. The resulting system of linear equations is solved with a direct method using the MUltifrontal Massively Parallel sparse direct Solver (MUMPS; Amestoy and others, 2001, 2006).

The current version of Elmer/Ice is not able to deal with a changing domain in the map plane. Thus the ice front is fixed in time, and a minimum ice thickness of 10m is applied everywhere and for all times. This implies that initially glaciated points are not allowed to become ice-free.

3. SeaRISE Experiments

3.1. Palaeoclimatic spin-up

In order to obtain a suitable present-day configuration of the Greenland ice sheet that can be used as an initial condition for future climate experiments, it is desirable to carry out a palaeoclimatic spin-up over at least a full glacial cycle. The resulting present-day conditions of the Greenland ice sheet can be particularly sensitive to the initialization method (Rogozhina and others, 2011). Here, similar to the spin-up described by Greve and others (2011), the forcing follows that specified by SeaRISE.

We have been unable to perform an entire spin-up with Elmer/Ice, due to the prohibitive computing time that would be required for such a long simulation. For this reason, we conduct the palaeoclimatic spin-up from 125 ka BP until 200 years BP with the shallow-ice model, SICOPOLIS. The horizontal resolution is 10 km and the vertical direction is discretized by 81 equidistant, terrain-following layers. After an initial relaxation over 100 years (starting from the present-day topography and isothermal conditions at−10°C everywhere) in order to avoid spurious noise in the computed velocity field (Calov, 1994), we keep the topography fixed over time in order to preserve a good fit between the simulated and observed present-day topographies.

The spin-up with SICOPOLIS is conducted only until 200 years BP because the initial conditions produced by a shallow-ice model then used in Elmer/Ice would produce an initial shock that influences the results obtained with the future climate experiments. In order to mitigate the effects of the initial conditions on the full Stokes model, two successive runs are conducted with Elmer/Ice to produce the present- day ice-sheet configuration. The first run is conducted from 200 to 100 years BP, starting with the SICOPOLIS output that is interpolated from the regularly spaced finite-difference grid of SICOPOLIS to the finite-element mesh of Elmer/Ice using a bilinear method. This run keeps the topography fixed and is intended to relax the initial shock originating from the switch from the shallow-ice to the full Stokes dynamics. The second run, from 100 years BP until the present, allows the ice-sheet surface to evolve, forced by a constant present- day climate. The obtained present-day configuration of the ice sheet is used as the initial condition for the future climate experiments with Elmer/Ice. By contrast, for the future climate experiments with SICOPOLIS, a fixed-topography spin-up with SICOPOLIS from 125 ka BP until the present is used.

3.2. Future climate experiments

For the future climate experiments, we use the same set of Sea RISE experiments as was employed by Greve and others (2011). This represents a subset of the suite defined in the ‘2011 Sensitivity Experiments’. Due to excessive computing times, we run them only for 100 rather than 500 years:

• Experiment C1: Constant climate control run; beginning at present (more precisely, the epoch 1 January 2004 0:00, corresponding to t = 0) and running for 100 years, holding the climate steady to the present climate.

• Experiment S1: Constant climate forcing with increased basal lubrication. This is implemented in Elmer/Ice by halving the basal drag (essentially doubling the basal sliding) everywhere in the domain.

• Experiment C2: AR4 climate run; starts with the same present-day condition, but the climatic forcing (mean annual temperature, mean July temperature, precipitation) is derived from an ensemble average from 18 of the AR4 models, run for the period 2004-98 under the A1B emission scenario; beyond 2098 the climate persists to the end of the run 100 years into the future. In order to avoid a sudden climate jump at the initial time, the forcing is applied by calculating precipitation and temperature anomalies (relative to 2004), which are then added to the present-day climate specified by Eqn (7) and the data by Ettema and others (2009).

• Experiment T1: Combination of C2 and S1, i.e. AR4 climate forcing with increased basal lubrication (halved basal drag).

In order to minimize the shock that arises from the transition from the fixed-topography spin-up to the future climate experiments with evolving topography, in neither case is the ice sheet allowed to extend beyond its present-day margin.

The remaining experiments defined in the ‘2011 Sensitivity Experiments’ of the SeaRISE group will be considered in future work.

4. Results and Discussion

The results of the initialization runs carried out with Elmer/Ice (Section 3.1) are shown in Figure 2. The surface topography is in good agreement with that observed (Bamber and others, 2001a, not shown), with differences of the order of tens of metres (see Table 2) as a consequence of the evolving free surface during the last 100 years. The surface velocity shows the expected distribution, with small velocities (<10ma^-1) around the major ice ridges and a general speed-up towards the coast. The pattern agrees well with Joughin and others’ (2010) interferometrically measured velocities, with the notable exception of the ‘Northeast Greenland Ice Stream’ (NEGIS) and the northwestern outlet glaciers where the velocities are relatively small. Basal temperatures are at pressure melting, most notably in the southwest and southeast, but also under the NEGIS. For the ice- core locations GRIP, NorthGRIP, Camp Century and Dye 3 where observations exist, Table 2 shows the comparison. The agreement is good for GRIP and Camp Century, but poor for NorthGRIP and Dye 3. This is probably due to shortcomings of the applied geothermal flux distribution by Shapiro and Ritzwoller (2004) and could be improved by the tuning method of Greve (2005); however, in this study we work with the Shapiro and Ritzwoller (2004) geothermal flux, according to the SeaRISE recommendation.

Fig. 2. Present-day configuration computed by Elmer/Ice starting from the SICOPOLIS palaeoclimatic fixed-topography spin-up at 200 years bp. (a) Surface topography, (b) surface velocity and (c) basal temperature relative to the pressure-melting point.

Table 2. Simulated (Elmer/Ice) and observed present-day ice thicknesses and basal temperatures for the ice-core locations GRIP, NorthGRIP, Camp Century and Dye 3.

Stable results could be obtained with both the full Stokes model, Elmer/Ice, and the shallow-ice model, SICOPOLIS, for all four future climate experiments described in Section 3.2. Simulated surface velocities for the control run, C1 (constant climate forcing), after 100 years are shown in Figure 3. The results for Elmer/Ice (Fig. 3a) show that major ice streams and outlet glaciers are active. In particular, Jakobs− havn Isbr? and Petermann, Kangerdlugssuaq, Helheim and further southeast outlet glaciers show continued fast flow. The NEGIS, however, is characterized by lower velocities of 30-100 m a^-1, with no pronounced acceleration towards the margins. The results for SICOPOLIS (Fig. 3b) exhibit generally higher ice-surface velocities around the ice margin. In addition to Jacobshavn Isbr? and Helheim and further southeast outlet glaciers, Kangerdlugssuaq Glacier, the outlet glaciers of the NEGIS and many further areas show fast flow with velocities exceeding 1000 m a^-1.

Fig. 3. Surface velocities (a, b) and basal temperatures relative to pressure melting (c, d) computed with Elmer/Ice (a, c) and SICOPOLIS (b, d) for experiment C1 (constant climate control run) at t = 100 years (year 2104).

This different dynamical behaviour of Elmer/Ice and SICOPOLIS near the ice margin has several causes. The representation of fast-flowing ice streams and outlet glaciers in Elmer/Ice benefits from the much finer grid resolution, and Elmer/Ice solves the full Stokes equations, so all components of the stress tensor are included. The consequence with respect to ice-stream dynamics is that Elmer/Ice accounts for the lateral drag resulting from local fast flow embedded in slower-flowing ice, which limits the velocity contrast, while the shallow-ice solver of SICOPOLIS does not exhibit lateral drag and thus tends to over-predict fast ice flow. Another reason for the generally lower surface velocities produced by Elmer/Ice lies in the different basal thermal conditions (Fig. 3c and d). The temperatures computed with Elmer/Ice after 100 years are generally lower and the temperate- based areas smaller, while the temperatures computed with SICOPOLIS are higher. The cooler conditions obtained with Elmer/Ice originate from the initial conditions (Fig. 2c) where the temperatures are generally lower than the initial conditions used by SICOPOLIS (not shown). However, the control run shows that the basal temperatures have generally increased in comparison to the initial conditions, particularly at Petermann and the northwestern outlet glaciers. This could indicate that the initial shock due to the sudden change of ice dynamics from shallow ice to full Stokes has gradually been smoothed out. Of course, the shallow-ice approximation used in SICOPOLIS applies also to Eqn (5) (neglect of the horizontal heat conduction), so the differences between the two models also result from the different temperature equations.

Figure 4 shows the results obtained for experiment S1 (doubled basal sliding). The surface velocities computed with both Elmer/Ice and SICOPOLIS show the expected response: an increase of the flow speed in all areas where the base is at or near the pressure-melting point. Consequently, both models produce faster-flowing ice streams and outlet glaciers compared to the control run, C1. The surface velocities computed with Elmer/Ice show higher sensitivities, with higher flow speeds observed at Jakobshavn Isbr? and the NEGIS and at the Petermann outlet glaciers. Elmer/Ice also produces more localized fast-flowing outlet glaciers at the northwestern margins. By contrast, the surface velocities computed with SICOPOLIS are only larger than their Elmer/Ice counterparts at the eastern margins, mainly at Kangerdlugssuaq Glacier, due to the larger area at the pressure-melting point. Here again the temperatures produced by Elmer/Ice are lower, but the increased basal heating, related to the larger basal sliding, allows the melting point to be reached at larger areas for the NEGIS and Jakobshavn Isbr? as well as the major outlet glaciers. It is also remarkable that although Elmer/Ice has smaller temperate- based areas than SICOPOLIS, the model shows a major speed-up of the ice-sheet flow, equal to or greater than that observed with SICOPOLIS. At the same time, for both models, the increased ice flow leads to increased advection of cold interior surface ice downwards and outwards, which should cool down the ice base compared to the control run. This is more evident for Elmer/Ice than for SICOPOLIS, perhaps due to the lower vertical resolution that does not capture so well the counteracting effect of increased strain heating near the base.

Fig. 4. Surface velocities (a, b) and basal temperatures relative to pressure melting (c, d) computed with Elmer/Ice (a, c) and SICOPOLIS (b, d) for experiment S1 (constant climate forcing, doubled basal sliding) at t = 100 years (year 2104).

The surface velocities and basal temperatures computed for experiments C2 (AR4 climate forcing) and T1 (AR4 climate forcing plus doubled basal sliding) are shown in Figures 5 and 6, respectively. For both Elmer/Ice and SICOPOLIS, the results are very similar to those obtained with experiments C1 and S1, respectively. Only a limited speed-up occurs because of the higher surface temperatures that penetrate slowly into the deeper ice. This indicates that the impact of a warmer climate in the absence of accompanying dynamical forcings has only a small effect on the dynamics and thermodynamics of the ice sheet on the considered timescale of 100 years, and acts mainly by the changed surface mass balance.

Fig. 5. Surface velocities (a, b) and basal temperatures relative to pressure melting (c, d) computed with Elmer/Ice (a, c) and SICOPOLIS (b, d) for experiment C2 (AR4 climate forcing) at t = 100 years (year 2104).

Fig. 6. Same as Figure 5, but for experiment T1 (AR4 climate forcing, doubled basal sliding).

Let us now focus on the detailed surface velocity evolution in the vicinity of Jakobshavn Isbr?. This is particularly interesting because the Greenland topography data used here (Section 2.2) are based on a special algorithm that preserves the continuity and depth of the trough below the ice stream (and below Helheim, Kangerdlussuaq and Petermann glaciers) in the gridded data (U.C. Herzfeld and others, unpublished information). Figures 7 and 8 show snapshots at t = 1, 10 and 100 years for experiments C2 and T1 conducted with Elmer/Ice and SICOPOLIS, respectively. On this zoomed spatial scale, the two sets of results are immediately distinguishable due to the coarser resolution employed by SICOPOLIS. Here the benefit of the mesh refinement manifests itself by a much smoother representation of the fast-flowing ice stream within the slower-flowing environment. In the Elmer/Ice results, an area of fast flow is visible and greatly expands out of the main bed trough at t =1 year for experiment T1. In the following (t =10 and 100 years), and for experiment C2, the velocities decrease and become more focused towards the margin. For experiment T1, the fast-flowing area outside the bed trough persists through time, with only a conspicuous and localized drop in velocity at t = 100 years. This local feature is due to a localized drop in the basal temperature which translates into a decrease in basal sliding. This sudden drop in temperature is probably not physical but rather related to some numerical issues. The SICOPOLIS results do not exhibit the widened areas of fast flow, because of the local nature of the shallow- ice approximation, whereas in Elmer/Ice the velocities are more influenced by the topography in the vicinity of the main trough, due to the non-local full Stokes force balance.

Fig. 7. Surface velocities in the area of Jakobshavn Isbræ computed with Elmer/Ice for experiments C2 (AR4 climate forcing) and T1 (AR4 climate forcing, doubled basal sliding) at t = 1 year (year 2005), 10 years (year 2014) and 100 years (year 2104).

Fig. 8. Same as Figure 7 but computed with SICOPOLIS.

The simulated ice volumes as functions of time are shown in Figure 9. The ice sheet reacts distinctly to all applied forcings. The control run, C1, with Elmer/Ice produces an ice volume gain of ~6cm s.l.e. during the 100 years of model time, while the same run with SICOPOLIS produces an ice volume loss of ~3 cm s.l.e. This difference, and in particular the stronger reaction of the Elmer/Ice run, is presumably due to the thermodynamically different initial conditions used by the two models. In order to largely remove this effect, we discuss the results of the three other experiments (S1, C2, T1) relative to the control run C1. After 100 years of model time, the ice volume losses, ΔV, are as follows:

Fig. 9. Ice volume (V) changes simulated with (a) Elmer/Ice and (b) SICOPOLIS for experiments C1 (constant climate control run), S1 (constant climate forcing, doubled basal sliding), C2 (AR4 climate forcing) and T1 (AR4 climate forcing, doubled basal sliding). Note that t = 0 corresponds to the year 2004.

The results from Elmer/Ice for the 2× basal sliding run, S1, show a ~43% higher sensitivity for ice volume loss than those from SICOPOLIS (computed as This is particularly remarkable because, as was discussed above, simulated flow velocities of Elmer/Ice are generally similar to SICOPOLIS, with only a few areas with faster velocities, so the higher sensitivity of Elmer/Ice is a consequence of the higher resolution of the fast-flowing areas. Moreover, Elmer/Ice seems to show that the full Stokes model has a higher sensitivity to dynamical changes, a crucial result for the implication of the Greenland destabilization due to increased basal lubrication.

For the direct global warming (AR4 climate) run, C2, Elmer/Ice is ~61% less sensitive than SICOPOLIS. The large difference is mainly explained by the initial conditions used by Elmer/Ice, where the present-day conditions computed by the spin-up run have higher surface elevations around the margins that built up during the last 100 spin-up years with freely evolving topography. This has an important consequence for the surface temperatures computed for the constant climate (Section 3.2) because the mean annual and mean July (summer) surface temperatures are dependent on the surface elevation (Eqn (7)), so the atmospheric lapse rate is effective. Because of the higher surface elevations at the margins for the initial conditions used by Elmer/Ice, surface temperatures are colder and surface melting is lower, implying that the evolution of the ice sheet from the imposed global warming driven by the change of the surface mass balance is more effective in the case of SICOPOLIS. The different sensitivities of Elmer/Ice and SICOPOLIS for experiment C2 are therefore not due to the differences between the two models (the representation of the surface mass balance is the same in both models), but mainly due to the different initial conditions.

For the combined (AR4 climate)/(2 × sliding) forcing of run T1, the sensitivities of both Elmer/Ice and SICOPOLIS are essentially equal to the sums of the sensitivities to the 2 × sliding forcing (S1) and the AR4 climate forcing (C2), which, in relative terms, makes Elmer/Ice ~21% more sensitive than SICOPOLIS. This near-linear behaviour results from the short model time of 100 years, during which the absolute changes in ice volume are limited to a few per cent, so the mutual influence between surface melting and ice flow remains small.

5. Conclusion

The full Stokes finite-element model Elmer/Ice was applied to the entire Greenland ice sheet. We carried out a set of SeaRISE experiments with it, and compared results with the SICOPOLIS shallow-ice model. This work marks an important step in ice-sheet modelling, as it is the first attempt to assess the likely range of the contribution of an ice sheet to sea-level rise in the future with a prognostic full Stokes model that captures ice dynamics most adequately.

Considering the computed surface velocities, the differences between the two force balances (full Stokes vs shallow ice) became evident. The surface velocities computed with Elmer/Ice are lower than those computed with SICOPOLIS for the control run and the AR4 climate run. For the fast flowing ice streams and outlet glaciers, this is mainly due to the lateral drag taken into account in the full Stokes model. The improved representation of fast-flow areas also benefited from the mesh refinement technique applied in Elmer/Ice that allows us to resolve them properly, while the regular 10 km grid of SICOPOLIS smears them out significantly. Further disparities were observed for the basal temperatures. The temperatures computed with Elmer/Ice are generally lower and the temperate-based areas smaller. This is possibly a shortcoming of the Elmer/Ice simulations that is due to the different initial conditions used by the models and the rather low vertical resolution of 17 layers (while SICOPOLIS employs 81 layers). So far we have not succeeded in increasing the vertical resolution because this leads to finite elements with a very small aspect ratio and thus numerical instabilities. For the experiments with dynamical forcing (S1 and T1), Elmer/Ice showed a higher sensitivity than SICOPOLIS with a greater acceleration of the ice flow at major ice streams and outlet glaciers. This greater speed-up in the case of Elmer/Ice produces surface velocities that are equal to or greater than the velocities obtained with SICOPOLIS.

The computed ice volume evolutions for the experiment with dynamical forcing (doubled basal sliding) showed a ~43% greater sensitivity for Elmer/Ice than for SICOPOLIS (relative to the constant climate control runs). The full Stokes approach of Elmer/Ice, along with the higher mesh resolution that leads to greatly improved representations of the fast-flowing zones, means that the model is more sensitive to dynamical destabilization processes, which is of great importance when investigating such phenomena to better estimate the resulting sea-level rise. Under the AR4 climate forcing, Elmer/Ice was ~61% less sensitive than SICOPOLIS, and under the combined forcing Elmer/Ice was ~21% more sensitive, in absolute terms essentially the sum of the two individual contributions. The higher sensitivity of SICOPOLIS for the AR4 climate forcing is mainly due to the different initial conditions used by the models; in the case of Elmer/Ice, the higher surface elevations near the ice-sheet margins limit surface melting.

Some important limitations of the results of this study must be noted. Although the full Stokes approach and the mesh refinement allow for an adequate representation of ice-stream dynamics, the applied Weertman-type sliding law, Eqn (8), is a severe simplification. It works reasonably well for the ice sheet as a whole (Greve and others, 1998; Greve, 2005), but its validity for fast-flowing ice with particular basal conditions is questionable. In addition, we have not attempted to account for the particular marine ice dynamics that play a role in changes of Jakobshavn Isbr?. This calls for improvements, to be tackled in future work. Most importantly, inverse methods will be applied in order to determine more suitable sliding laws for the ice streams and outlet glaciers, to be constrained by interferometrically measured present-day surface velocities (Joughin and others, 2010). This will also require improving the initial conditions and running Elmer/Ice in full Stokes for a larger part, or even the duration, of the spin-up. Further desirable changes concern treatment of the ice margins (implementation of moving margins, temperature boundary conditions) and a higher vertical resolution (which requires overcoming numerical stability issues due to extremely shallow finite elements). These improvements are partly underway and will hopefully lead to more accurate and reliable simulations of ice volume variations under changing climates.