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 (Reference Lemke and SolomonLemke and others, 2007; Reference Rignot, Velicogna, Van den Broeke, Monaghan and LenaertsRignot 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 (Reference Zwally, Abdalati, Herring, Larson, Saba and SteffenZwally and others, 2002; Reference Rignot and KanagaratnamRignot and Kanagaratnam, 2006; Reference Howat, Joughin and ScambosHowat and others, 2007; Reference JoughinJoughin 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 (Reference HutterHutter, 1983; Reference MorlandMorland, 1984) has been widely used for ice-sheet models (e.g. Reference HuybrechtsHuybrechts, 1990; Reference Calov and HutterCalov and Hutter, 1996; Reference Ritz, Rommelaere and DumasRitz and others, 2001; Reference Saito and Abe-OuchiSaito and Abe-Ouchi, 2004; Reference Rutt, Hagdorn, Hulton and PayneRutt 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. Reference Greve and BlatterGreve and Blatter, 2009). This gave rise to the formulation of higher-order models (Reference BlatterBlatter, 1995; Reference Baral, Hutter and GreveBaral and others, 2001; Reference HindmarshHindmarsh, 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 (Reference Dahl-JensenDahl-Jensen, 1989; Reference BlatterBlatter, 1995; Reference Colinge and BlatterColinge and Blatter, 1998; Reference PattynPattyn, 2000), and Reference Dahl-JensenDahl-Jensen (1989) demonstrated the importance of longitudinal deviatoric stresses for plane flow along a flowline. Reference PattynPattyn (1996, Reference Pattyn2000) and Reference Pattyn and DecleirPattyn and Decleir (1998) applied a 2-D higher-order model with thermomechanical coupling to Shirase drainage basin, Dronning Maud Land, Antarctica. Reference PattynPattyn (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 (Reference Huybrechts and PayneHuybrechts and others, 1996; Reference PaynePayne and others, 2000) and provided a comparison with the Saito-Blatter model that also includes higher-order dynamics (Reference Saito, Abe-Ouchi and BlatterSaito and others, 2003). More recent developments are the models of Reference Pollard, DeConto, Hambrey, Christoffersen, Glasser and HubbardPollard and DeConto (2007, Reference Pollard and DeConto2009) and Reference Bueler and BrownBueler and Brown (2009) that employ heuristic combinations of the shallow-ice and shallow-shelf approximations (Reference Morland, Van der Veen and OerlemansMorland, 1987; Reference MacAyealMacAyeal, 1989), as well as the application of a first- order model to the Greenland ice sheet by Reference Price, Payne, Howat and SmithPrice 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 (Reference SolomonSolomon 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. Reference GudmundssonGudmundsson, 1999; Reference Sugiyama, Gudmundsson and HelbingSugiyama and others, 2003; Reference Martín, Navarro, Otero, Cuadrado and CorcueraMartfn and others, 2004; Reference Price, Waddington and ConwayPrice and others, 2007; Reference Jouvet, Huss, Blatter, Picasso and RappazJouvet 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) (Reference PattynPattyn 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 Reference Ren and LeslieRen and Leslie (2011) and Reference Ren, Fu, Leslie, Chen, Wilson and KarolyRen and others (2011a, Reference Ren, Fu, Leslie, Chen, Wilson and Karolyb), 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. Reference Zwinger, Greve, Gagliardini, Shiraiwa and LylyZwinger and others, 2007; Reference Gagliardini and ZwingerGagliardini and Zwinger, 2008; Reference Durand, Gagliardini, Zwinger, Le Meur and Hind-marshDurand and others, 2009; Reference Zwinger and MooreZwinger and Moore, 2009; Reference Seddik, Greve, Zwinger and PlacidiSeddik 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 Reference Greve, Saito and Abe-OuchiGreve 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. Reference GreveGreve, 1997, Reference Greve, Weis and Hutter2000; Reference Greve, Saito and Abe-OuchiGreve 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. Reference Greve and BlatterGreve 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. Reference Greve and BlatterGreve 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 Reference Greve, Saito and Abe-OuchiGreve 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 Reference Fausto, Ahlstrom, Van As, Boggild and JohnsenFausto 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 (Reference DansgaardDansgaard and others, 1993; Reference JohnsenJohnsen 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 Reference EttemaEttema 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 Reference ReehReeh’s (1991) positive degree-day (PDD) method, supplemented by the semianalytical solution for the PDD integral by Reference Calov and GreveCalov 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 (Reference Huybrechts and de WoldeHuybrechts and de Wolde, 1999). Furthermore, the standard deviation of short-term, statistical air-temperature fluctuations is σ = 5°C (Reference Huybrechts and de WoldeHuybrechts and de Wolde, 1999), and the saturation factor for the formation of superimposed ice is chosen as P _max = 0.6 (Reference ReehReeh, 1991). Conversion from the present-day mean annual precipitation (Reference EttemaEttema and others, 2009) to the snowfall rate (solid precipitation) is done on a monthly basis using the empirical relation of Reference MarsiatMarsiat (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 (Reference GreveGreve, 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 (Reference Hindmarsh and Le MeurHindmarsh 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. Reference Greve and BlatterGreve and Blatter, 2009), where the geothermal flux distribution over the bedrock is given by Reference Shapiro and RitzwollerShapiro 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 Reference Joughin, Smith, Howat, Scambos and MoonJoughin and others (2010) with gaps filled by balance velocities of Reference Bamber, Layberry and GogineniBamber 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 Reference Frey and AlauzetFrey and Alauzet (2005), and the adaptation is carried out with the automatic tool YAMS (Reference FreyFrey, 2001; Reference Morlighem, Rignot, Seroussi, Larour, Ben Dhia and AubryMorlighem 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 Reference Joughin, Smith, Howat, Scambos and MoonJoughin and others (2010) with gaps filled by balance velocities of Reference Bamber, Layberry and GogineniBamber 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 (Reference Franca and FreyFranca and Frey, 1992; Reference Franca and FreyFranca 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; Reference Amestoy, Duff, L’Excellent and KosterAmestoy and others, 2001, Reference Amestoy, Guermouche, L’Excellent and Pralet2006).

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

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 (Reference Bamber, Ekholm and KrabillBamber 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 Reference Joughin, Smith, Howat, Scambos and MoonJoughin 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 Reference Shapiro and RitzwollerShapiro and Ritzwoller (2004) and could be improved by the tuning method of Reference GreveGreve (2005); however, in this study we work with the Reference Shapiro and RitzwollerShapiro 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 (Reference Greve, Weis and HutterGreve and others, 1998; Reference GreveGreve, 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 (Reference Joughin, Smith, Howat, Scambos and MoonJoughin 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.