Ice-dynamic projections of the Greenland ice sheet in response to atmospheric and oceanic warming

Continuing global warming will have a strong impact on the Greenland ice sheet in the coming centuries. During the last decade (2000–2010), both increased melt-water runoff and enhanced ice discharge from calving glaciers have contributed 0.6±0.1 mm yr to global sea-level rise, with a relative contribution of 60 and 40 % respectively. Here we use a higher-order ice flow model, spun up to present day, to simulate future ice volume changes driven by both atmospheric and oceanic temperature changes. For these projections, the flow model accounts for runoff-induced basal lubrication and ocean warming-induced discharge increase at the marine margins. For a suite of 10 atmosphere and ocean general circulation models and four representative concentration pathway scenarios, the projected sea-level rise between 2000 and 2100 lies in the range of+1.4 to+16.6 cm. For two low emission scenarios, the projections are conducted up to 2300. Ice loss rates are found to abate for the most favourable scenario where the warming peaks in this century, allowing the ice sheet to maintain a geometry close to the presentday state. For the other moderate scenario, loss rates remain at a constant level over 300 years. In any scenario, volume loss is predominantly caused by increased surface melting as the contribution from enhanced ice discharge decreases over time and is self-limited by thinning and retreat of the marine margin, reducing the ice–ocean contact area. As confirmed by other studies, we find that the effect of enhanced basal lubrication on the volume evolution is negligible on centennial timescales. Our projections show that the observed rates of volume change over the last decades cannot simply be extrapolated over the 21st century on account of a different balance of processes causing ice loss over time. Our results also indicate that the largest source of uncertainty arises from the surface mass balance and the underlying climate change projections, not from ice dynamics.


Introduction
Volume changes of the Greenland ice sheet result from a balance between ice accumulation on its surface and ice loss around its margin by both melt-water runoff and ice discharge into the surrounding ocean. In the thirty year period prior to 1990, the ice sheet geometry has been in a virtual balance with the prevailing climate but has 5 since been losing mass at an increasing rate (Rignot et al., 2011;Zwally et al., 2011;Shepherd et al., 2012;Sasgen et al., 2012). Almost half of the recent mass loss is attributed to increased ice discharge at the marine margins (van den Broeke et al., 2009;Shepherd et al., 2012;Sasgen et al., 2012;Vaughan et al., 2013). In the period 1972 to 1995, glacier terminus positions and ice flow were rather stable around Greenland 10 (Howat and Eddy, 2011;Bevan et al., 2012). Over the last decade, however, ice sheetwide surface velocity observations reveal complex spatial and temporal changes with accelerated glacier flow in the northwest, more variability in the southeast and relatively steady flow elsewhere (Moon et al., 2012). Since 2000, glaciers in the southeast showed synchronous retreat and acceleration followed by a general deceleration and 15 re-advance after 2005 (Howat et al., 2007;Moon and Joughin, 2008;Murray et al., 2010). In the northwest, the dynamic glacier response was more asynchronous but most outlet glaciers have accelerated and retreated at some point between and 2009(McFadden et al., 2011Moon et al., 2012).
A prominent example for recent dynamic changes of outlet glaciers in west Green- 20 land is Jakobshavn Isbrae. Starting in 1998, its frontal zone sped up from about 6 to 12 km yr −1 within 5 years (Joughin et al., 2004(Joughin et al., , 2008c. This speedup coincided with an upstream propagation of glacier thinning by increased flow divergence. One hypothesis links the main acceleration to a successive loss of buttressing on the grounded ice as the floating ice tongue destabilised and collapsed while another hypothesis points 25 to a speedup initiated by a weakening of the ice in the lateral glacier margins (van der Veen et al., 2011). In any case, the initiation of the glacier acceleration and retreat co- incides with an intrusion of warm Atlantic Water into Disco Bay that likely entered the local fjord systems (Holland et al., 2008). In southeast Greenland, speedup and retreat peaked in 2005 for Helheim and Kangerlussuaq Glaciers, which are both located at the end of ∼ 80 km long fjords. Though the acceleration peak occurred simultaneously for both glaciers, the speed and 5 retreat behaviour leading to this event differs (Stearns and Hamilton, 2007;Joughin et al., 2008b). While Helheim shows a continuous acceleration starting in 2002 with a cumulative retreat of the ice front of 8 km, Kangerlussuaq exhibited an abrupt retreat and acceleration between 2004 and 2005. Yet for both glaciers, the acceleration events were temporary and glacier speeds dropped again to the pre-speedup level 10 (Bevan et al., 2012). In recent years, there is growing evidence that changes in fjord circulation and fjord stratification were the first-order control on this regional retreat and acceleration (Murray et al., 2010;Straneo et al., 2011). The southeast glacier acceleration was preceded by a period of low runoff from the ice sheet, which weakened the East Greenland Coastal Current and allowed warm waters from the Irminger Sea 15 to reach the coast. Subsequently, the coastal current regained strength and provided again cold arctic waters, which presumably has led to the regional re-stabilisation.
At the northern margin of the Greenland ice sheet, the Petermann glacier recently lost a major part of its 80 km long floating tongue. On 4 August 2010, about one fifth of the ice tongue broke off and drifted out of the fjord into Nares Strait (Falkner et al., 20 2011). In line with the above speedup examples, this breakup event was also preceded by ocean warming in the hundred meters above the 300 m deep sill at the southern end of Nares Strait (Münchow et al., 2011).
Warm and saline waters with tropical origin are in fact found at intermediate depth beyond the continental shelf break all around Greenland. There is evidence that these 25 waters can overcome the sills at the mouth of local fjord systems around Greenland (Straneo et al., 2012). Warming of deep fjord water can intensify submarine melt below an existing ice shelf or mélange cover (Motyka et al., 2011), or directly at the calving front (Rignot et al., 2010). Thinning in the frontal zone can, in turn, reduce the buttress-3854 base (Schoof, 2010). Observations on both ice velocity and local runoff at various positions along the western flank of the Greenland ice sheet show a speedup in summer when surface meltwater production peaks (Zwally et al., 2002;van de Wal et al., 2008;Bartholomew et al., 2011;Sundal et al., 2011;Shannon et al., 2013). The total amount of melt-water production then determines the local effect on the annual acceleration 10 above the winter reference (Sundal et al., 2011). Basal lubrication is also hypothesised to enhance ice flow towards the marine margin and to thereby influence ice discharge.
While ice discharge changes explain about 40 % of the recent ice loss on Greenland, the remainder is attributed to a decreasing surface mass balance (SMB; Vaughan et al., 2013). Most direct observations of the SMB components have locat at most regional 15 character and are limited to the last decade (van den Broeke et al., 2011). Therefore, they are too short and not representative to directly infer ice-sheet wide trends. Yet SMB modelling has improved with the availability of validation data. Regional climate models are now capable of producing a physically-based ice-sheet wide SMB estimate (Vernon et al., 2013). SMB model results show that the five years with highest annual 20 meltwater runoff since 1870 fall into the period after 1998 (Hanna et al., 2011). This concentrated occurrence of years with peak runoff indicates the general increase in runoff over the last decades, which explains the negative trend in the SMB since the late 1990s (Ettema et al., 2009). In addition, the cumulated melt area has continuously increased, and melt extents since 2000 are on average twice as large as in the early 25 1980s (Fettweis et al., 2011. For ice loss on Greenland over the next few centuries, a major contribution is expected from a decreasing SMB, or more precisely an increase in surface meltwater runoff (Cazenave et al., 2013). By now, the modelling community has achieved to im- Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | prove regional climate models (RCM) to the point that they reproduce past and present changes in various components of the SMB (Vernon et al., 2013). Owing to a shortage in the observational coverage, the largest source of model uncertainties remains in the treatment and quantification of meltwater percolation and refreezing within the snowpack. The physical complexity of RCMs often limit any application on ice sheet-wide 5 scales to the usage of coarse resolutions (beyond 10 km). Yet it is within a narrow band of several tens of kilometres around the ice sheet margin that the largest SMB changes are expected under atmospheric warming. Assuming small perturbations, RCM simulations often use a fixed ice-sheet geometry. But under strong future warming, margin thinning attains a level, where SMB models need to account for these geometric 10 changes. Instead of using a downsampling procedure for RCM SMB fields (Franco et al., 2012), parametric SMB approaches are often applied in high-resolution ice-flow models (Huybrechts, 2002;Robinson et al., 2011;Greve et al., 2011). Though such approaches rely on parameterisations of individual SMB components, volume projections can account for the feedback between geometric adjustments and SMB changes. 15 This feedback comprises the effect of changes in ice discharge on the SMB.
Here we make an effort towards including more ice-dynamical processes in a thermomechanically-coupled, three-dimensional ice flow model (Huybrechts and de Wolde, 1999) with the aim to better assess the impact of ice dynamics on ice volume projections. These projections are driven by the Representative Concentration Pathways 20 (RCP; Moss et al., 2010) used for the IPCC's Fifth Assessment Report (AR5) (refer to IPCC, 2013). The ice dynamic model component includes parameterisations for ocean warming-induced discharge increase and runoff-induced basal lubrication (Sect. 2). A selection of ten Atmosphere and Ocean General Circulation Models (AOGCM) in the CMIP5 dataset (Taylor et al., 2012) provides the climatic input for the projections. 25 From this input, both atmospheric and oceanic forcing is prepared as anomalies to drive the ice sheet model (Sect. 3). In order to be able to cover the full range of scenarios and climate models, a parametric SMB model is applied. Its ablation component is based on the positive degree-day (PDD) method that accounts for retention and re-  (Janssens and Huybrechts, 2000). Before the actual sea-level projections of the Greenland ice sheet are presented in Sect. 5, the model evolution in the recent past is used for validation against direct and indirect observations (Sect. 4).

Model description and initialisation
The three-dimensional thermo-mechanically coupled ice sheet model comprises three 5 main components that respectively describe the mass balance at the upper and lower ice sheet boundaries, the ice dynamic behaviour and the isostatic adjustment of the Earth lithosphere (Huybrechts and de Wolde, 1999;Huybrechts, 2002;Fürst et al., 2011).
The simulated ice flow arises as a viscous response of the material to gravitational 10 forcing. Using a higher-order approximation to the force balance, the model accounts for effects from horizontal gradients in membrane stresses (Fürst et al., 2011). This ice-dynamic core allows for a more realistic inland transmission of perturbations at the ice sheet margin . The flow component of the ice sheet model also accounts for the direct effect of ocean warming on ice discharge and for runoff-induced 15 lubrication. Both effects are parameterised and presented in the following sections.
In the parametric SMB model, a distinction is made between snow accumulation, meltwater runoff and meltwater retention in the snowpack. Net surface accumulation is based on the Bales et al. (2009) map for the present ice accumulation representative for the period 1950-2000. For the ablation component, the melt and runoff model re-20 lies on the widely used positive degree-day runoff/retention approach (Janssens and Huybrechts, 2000;Gregory and Huybrechts, 2006). This approach first determines the positive degree-day sum from monthly air temperature input assuming a statistical variability in daily near-surface temperatures around the monthly mean (with a standard deviation of 4.2 • C). Melt rates are then determined with degree-day factors for snow 25 and ice melting of respectively 0.0030 and 0.0079 m d −1 • C −1 ice equivalent (i.e.). Initial melt in snow-covered regions is stored as capillary water until the snowpack becomes ice sheet model grid. The geothermal heat flux is inferred from seismic data (Shapiro and Ritzwoller, 2004) and was adjusted with Gaussian functions at the deep ice core sites (NEEM, GRIP, NGRIP, Dye3 and Camp Century) to reproduce observed basal temperatures. 20 Observations on ice velocities show seasonal speed-up in the summer period when surface melt-water production peaks (Zwally et al., 2002;Rignot et al., 2010;Bartholomew et al., 2011). Surface runoff generally finds a way into the ice body through moulins and the water is assumed to reach the bed near the ice-sheet margin. The rate of basal melt-water discharge determines the two-fold character of the 25 subglacial drainage system, which in turn controls lubrication and its effect on the sliding velocity (Schoof, 2010). Though this mechanism has a seasonal cycle, it is the integrated effect over one year and changes therein which are important for the future 3858 Introduction ice-dynamic response. Therefore, we assume a relation between the annual surface runoff and the annual, relative sliding increase with respect to the winter reference. In the ice flow model, the Weertman sliding relation is therefore extended with a multiplier S BL that depends on the annual rate of basal meltwater discharge.

Effect of surface runoff on basal lubrication
Here sliding velocities are denoted with v b , basal drag with τ b , the sliding factor with A S and the ice thickness H. In this parameterisation, the basal meltwater discharge rate is assumed equal to the local surface runoff R, whilst neglecting contributions from basal melting or meltwater routing beneath the ice sheet. Theoretical work on 10 subglacial drainage systems indicates a speedup peak for a specific rate of basal water discharge (Schoof, 2010). Above this discharge rate, a channelised basal drainage system is preferred, which is associated with lower relative speed-up values. In the absence of local runoff, no lubrication effect is simulated. These characteristics are approximated by a Poisson-like functional dependence between relative speedup and 15 runoff. (2) In this notation, the unknown parameters a, b and c are assumed positive. Instead of assessing the uncertainties associated with the exact parameter choice in this func-20 tional dependence, as in Shannon et al. (2013), one set of parameters is determined and applied for the projections.
The three unknown parameters are determined using observational data on annual acceleration and runoff at two locations along the western flank of the Greenland ice sheet (Fig. 1). The first location is east of Kangerlussuaq and upstream of Russell 25 glacier often referred to as the K-transect (van de Wal et al., 2008;Shepherd et al., 2009;Bartholomew et al., 2010Bartholomew et al., , 2011. Here a consistent picture emerges with annual accelerations of up to 20 % above the winter background for runoff rates below 3.5 m 3859 Introduction observed and modelled critical runoff rates needs to be considered in our functional dependence ( Fig. 1). For a modelled ice sheet geometry, the used mass balance model gives annual runoff rates of up to 4 m yr −1 near the K-transect. But due to a faster inland decrease in modelled runoff, as compared to observations, upstream speedup would be underestimated. Taking this into account, the following parameter values are cho- about 1 m yr −1 . At these locations however, the annual speedup is also influenced by seasonal changes at the marine termini. These observational estimates are covered by the suggested functional dependence. The presented parameterisation might be affected by the observational bias towards the western flank of the Greenland ice sheet. Yet the approach accounts for difference between observed runoff estimates and sim- 25 ulated runoff values from a model that is calibrated to reproduce the ice sheet-wide SMB.

Effect of ocean warming on ice discharge
With the aim to parameterise ocean-induced changes in ice discharge, outlet glacier accelerations are linked to oceanic warming assuming a uniform functional dependence. This relation is derived by relating interferometric velocity observations (Rignot and Kanagaratnam, 2006;Moon et al., 2012) to temperature variability diagnosed from 5 five ocean basins in available AOGCMs for the decade 2000-2010. Observations during this decade show an average speedup of outlet glaciers in the southeast of 34 % and in the northwest of 28 %, while other regions show no significant trend (Moon et al., 2012). Scaling these accelerations to the entire ice sheet and weighting them with the regional discharge distribution (Rignot and Kanagaratnam, 2006) results in an aver-10 age ice discharge increase of about 10 to 15 %. This increase shows an almost linear trend over the last decade (Rignot et al., 2011). Accounting for additional discharge estimates inferred from observations and regional climate models (Sasgen et al., 2012), the decadal increase in ice discharge explains between 25 and 40 % of the total mass loss (Shepherd et al., 2012). Considering the oceanic temperature forcing at hand to- 15 gether with the fast marginal adjustment properties of the ice sheet model , a linear increase in discharge is best simulated by a non-linear relation between ocean temperatures and sliding velocities. In addition, results from a generalisation of the flow-line response of individual outlet glaciers to a large-scale Greenland ice-sheet application Goelzer et al., 2013) support the choice for a exponen-20 tial dependence. The selected relationship is calibrated such that the ice sheet model reproduces the relative contribution of the discharge increase to the total ice loss over the last decade in response to the considered climate models.
The amplification of the sliding factor A outlet S applies exclusively to marine-terminated glaciers using the temperature anomaly ∆T ocean in the adjacent ocean basins. The forcing is applied up to 20 km inland from the calving front for ice grounded below sealevel to account for far-reaching loss in back-stress (Nick et al., 2012).
Prescribing an experiment with a linear increase in ocean temperatures under constant atmospheric forcing, the ice sheet model shows an increase in ice discharge that is more regular than the amplification of the sliding coefficient (Fig. 2). The reason is 5 a geometric adjustment and thinning at the marine margins that limits the attainable ice export. After 100 years, ocean temperatures are kept constant and ice discharge remains at an increased level. Yet the ongoing geometric adjustment causes a general decrease of the ice discharge in this latter period. 10 As initialisation to the present-day state, the model is spun up over a full glacial cycle as described in Huybrechts (2002). The ice sheet geometry evolves freely in response to past changes in regional surface temperatures, in precipitation and in sea-level stand. Though the general approach is unchanged, the underlying reconstruction for past temperature changes is updated with recent proxy information from several ice cores 15 (for details see Appendix A). A new compilation of accumulation observations over the Greenland ice sheet (Bales et al., 2009) is used as basis for scaling past precipitation changes with the mean annual temperature change (by 5 % • C −1 ). Finally, a new parameterisation to improve the retreat history from the Last Glacial Maximum is applied (Simpson et al., 2009), which is constrained by proxies for relative sea level. Experi-20 ments have shown that switching at 3 kyr BP from a shallow ice approximation to the higher-order formulation is sufficient to resolve the main effects of including horizontal stress gradients. Using an unconstrained model evolution during the spin-up phase guarantees a selfconsistent model state in the present day but the geometry will deviate from the 25 observed state. Therefore, key model parameters are tuned to minimise geometric and dynamic differences after the spin-up. To efficiently sample the parameter space, a Latin hypercube sampling (LHS; McKay et al., 1979)  nique was previously used for assessing the parameter sensitivity of an ice sheet model spin-up (Stone et al., 2010). But instead of only varying the positive degree-day factors both for ice (DDF ice ) and snow (DDF snow ) together with an enhancement factor (m) to the rate factor, one additional parameter is included in the sampling here: i.e. the sliding coefficient (A S ). These four parameters control both the SMB and the dynamic 5 state of the modelled ice sheet. Parameters are selected in a range of 75-125 % for the degree-day factors, 36-450 % for the enhancement factor m and 50-200 % for A S with respect to a previous calibration (Table 1). Eight criteria were chosen to quantify differences between the modelled ice sheet and the observed present-day state. The minimisation includes differences in the total 10 ice volume, in the ice-covered area, in the ice area above 3000 m and below 1500 m surface elevation, in the southwest position of the land-terminated ice margin and in the global ice thickness and surface elevation. Instead of exclusively focussing on geometric tuning diagnostics, as in (Stone et al., 2010), a final criterion evaluates the dynamic state of the ice sheet. Ice discharge in the decades prior to 1990 is assumed to have 15 compensated for ∼ 60 % of the average accumulation (Ettema et al., 2009). This additional criterion considerably reduces the parameter space. But all criteria were considered equally when assessing each parameter combination. The parameter set with the overall best performance (Table 1) shows very similar positive degree-day factors than the previous tuning but parameters for ice flow are slightly reduced. This reduction is 20 necessary because of higher velocities in the ablation zone when using the parameterisation for runoff-induced speedup. Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | meteorological reanalysis and ECMWF operational analysis data as described in Hanna et al. (2011). Anomalies and ratios are calculated with respect to the period . This assumes that the ice sheet was in quasi-equilibrium with the prevailing climate of that time, as is often done (e.g. Hanna et al., 2005). The reference precipitation distribution is from (Bales et al., 2009). In the same way, the oceanic temperature 5 anomalies are calculated from the Atmosphere and Ocean General Circulation Models (AOGCMs).

Future scenarios
For future ice sheet simulations, climate projection data from ten AOGCMs were selected from the WCRP's CMIP5 multi-model dataset prepared for the IPCC AR5 (Taylor  (Table B1 gives a complete overview of the considered AOGCMs). For these projections, the AOGCMs were forced with four CMIP5 Representative Concentration Pathway (RCP) scenarios (Moss et al., 2010). 15 The same anomaly approach as for the reference period is used to avoid a bias by the mean states of the AOGCMs. Monthly surface air temperature anomalies, annual precipitation ratios and annual ocean temperature anomalies are therefore considered with respect to the same 1960-1990 reference period. 20 Monthly surface air temperature anomalies and annual precipitations ratios are derived for each individual AOGCM over the ice sheet model domain. These future atmospheric anomalies drive the SMB model starting from the year 2005. In most cases, the data covers the period up to 2100 or 2300 AD. Missing data in the last year of two AOGCMs were filled by repeating the previous year. The annual air temperature anomaly averaged over the present ice-sheet extent ( Fig. 3) is instructive as a general trend but conceals the 2-D pattern of the warming (not shown). In general, the spatial pattern of the temperature forcing shows an expressed north-south gradient of up to 10 • C by 2100 AD. This latitudinal gradient depends on the climate sensitivity and the polar amplification of each AOGCM. For one 5 latitude, a more expressed warming on either side of Greenland depends highly on the individual AOGCM. In general, the precipitation increase over the model domain scales with the scenario intensity. By 2100, the ensemble averages per RCP show 13, 19, 23 and 37 % additional precipitation for respectively RCP2.6, RCP4.5, RCP6.0 and RCP8.5. For RCP2.6 and RCP4.5, these values increase to respectively 19 and 31 % 10 by 2300.

Ocean forcing
Oceanic forcing is decomposed into time series for five different oceanic basins. projections, ocean temperature changes in these basins are related to ice discharge changes at the marine-terminated margin of the Greenland ice sheet. Ocean circulation in the deeper ocean around Greenland, off the continental shelf, is resolved in most AOGCMs. Ocean basins are latitudinally delineated by the 60 • N,

70
• N, 80 • N parallels and the Northpole at 90 • N, and confined by the Greenland coast-5 line (Fig. 4). In each individual basin, AOGCM grid box centres that lie within a 300 km radius from the Greenland coastline are considered. This belt covers the continental shelf and a part of the deep ocean beyond the shelf break. The resulting basin temperature anomalies are not very sensitive to a radius increase to 500 km. In the vertical, temperatures are averaged over a depth of 200 to 600 m. The upper limit is inspired 10 by the average freshwater layer thickness in Greenlandic fjords (Straneo et al., 2010(Straneo et al., , 2011(Straneo et al., , 2012 together with intermediate depth locations of offshore AW (Holland et al., 2008). The latter argument combined with the fact that Greenlandic fjords have typical sill depths of several hundred metres gives rise to the lower bound. Area-and depthaveraging of all AOGCM grid points in each basin provides five temperature time series 15 for each AOGCM and each RCP. Ocean temperature anomalies for each basin are considered with respect to the 1960-1990 average (Fig. 5). For each basin, the annual temperature anomaly records are filtered with a 5 year moving average. This is necessary to prevent high frequency oscillations when forcing the ice-dynamic model. Though there is a tendency 20 for stronger warming in the northern ocean basins in many of the AOGCMs, differing trends within the five basins are highly dependent on the individual climate model.

Ice sheet evolution in the recent past
After the glacial-cycle spin-up, the present-day ice sheet geometry is in a selfconsistent state concerning ice geometry, dynamics, temperature and SMB. The geom-25 etry and temperature naturally carry the long-term memory of the ice sheet evolution. The main shortcoming from such a spin-up is that for the present day the modelled ge- Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | ometry does not exactly match observations. Like in other studies with a similar spin-up technique, ice thicknesses near the margin tend to be over overestimated and therefore the ice extent is somewhat larger (e.g. Huybrechts, 2002;Robinson et al., 2011;Greve et al., 2011;Graversen et al., 2011). Though the geometric mismatch biases the SMB near the margin, the ice sheet-wide SMB compares well with other approaches 5 (see below). Thicker margins also affect the modelled ice flow, which then tends to underestimate margin ice velocities (Fig. 6). A side-by-side comparison shows that the locations and the magnitudes of channelised ice flow towards the marine margin are fairly well reproduced on the 5 km grid. In this spin-up technique, regions of fast flow naturally arise from the interplay between deformation, sliding and thermo-dynamics. 10 Since velocities generally drop below 100 m yr −1 within some tens of kilometres from the ice-sheet margin, regions further upstream are not expected to directly contribute to ice discharge within one century. More meaningful than matching velocities at the margin is therefore that the model is capable of reproducing ice discharge rates and their regional distribution around Greenland ( Table 2). The simulated present-day state 15 shows a total ice discharge that slightly exceeds otherwise inferred values. The 5 % overestimation mostly arises from simulated ice-ocean contact in regions where no ice-sheet cover is observed, i.e. in the north and east. On 5 km resolution, ice flow towards the margin is more channelised and the regional agreement between modelled and inferred discharge improves, on a regional level and down to the level of major out-20 let glaciers. The match on a drainage basin level arises naturally without any additional model tuning. A 20 km model spin-up is only capable of reproducing the large-scale regional distribution and the total ice discharge. In this regard, the glacial-cycle spin-up method is preferable to another initialisation technique that aims at inverting for observed ice velocities using the observed geometry (Gillet-Chaulet et al., 2012). Though 25 it reproduces observed velocities, this latter initialisation technique is confronted with a strong initial model drift. Therefore, we believe that the free-geometry spin-up, using a model with increased dynamic complexity on high resolution, provides a useful ble produces the inferred ∼ 40 % share of the total mass loss. However, the mean is at the lower end of observations during this period and results from a weak oceanic warming around Greenland over the last decade in the used climate models (Fig. 3). Figure 7 and Table 4 summarise the volume projections of the Greenland ice sheet 10 for all models and all scenarios under investigation. A breakdown by individual climate models is presented in Appendix B. By 2100 AD, the full model and scenario range of Greenland sea-level contributions is between 1.4 and 16.6 cm ( Fig. 7 and Table B1). This range is slightly higher than the 1-12 cm found for the IPCC AR4 (Meehl et al., 2007), which included the additional uncertainty arising from the SMB model.

15
The higher maximum in sea-level projections is somewhat unexpected because the RCP scenarios have a reduced upper bound for radiative forcing by 2100, when compared to the previously used scenarios. Yet the larger range is attributed to accounting for the future discharge increase in these simulations. In terms of the SMB contribution to future ice loss, the IPCC AR5 (Cazenave et al., 2013) gives a range of 1-11 cm, 20 confirming the results of the previous report. But the AR5 comprises first attempts to quantify the contribution from future changes in ice discharge. It states an additional contribution from dynamic changes of 1-9 cm for all RCP scenarios. The new AR5 suffers, however, from a multitude of studies with not directly comparable forcing or setups when quantifying feedbacks between ice dynamics and surface mass balance. Their Until 2050 AD, there is hardly any difference in the mean sea-level contribution between the four scenarios. This is in agreement with a similar behaviour for the underlying atmospheric and oceanic forcing (Sect. 3.2.1). The ensemble spread in sea-level evolution for each scenario arises from the different climate trajectories followed by the individual AOGCMs. This spread is largely overlapping during the first century for three 5 scenarios. The exception is RCP8.5, a high-impact scenario assuming a high-emission fossil-fuel orientated world. This scenario causes a mean centennial sea-level contribution of 10.2 cm, which is about twice as large as for the other RCPs. The reason is an average warming of ∼ 7 • C over Greenland that is also more than twice as high. In addition, RCP8.5 is the only scenario for which mass loss rates significantly increase 10 throughout the next century.
Projections for the two lowest scenarios were continued until 2300 AD. Both assume a stringent climate policy with focus either on terrestrial carbon for mitigation (RCP4.5) or on negative emissions (RCP2.6). Both scenarios aim for a climate stabilisation but only RCP2.6 has a peak greenhouse gas concentration before 2100 AD and declines 15 afterwards (Moss et al., 2010). For both scenarios, the Greenland contribution to global sea-level rise increases continuously but for RCP2.6 the rate of increase gradually levels off. In this case, it appears that a new ice-sheet equilibrium with limited ice loss (< 20 cm of sea-level rise) is attainable. For RCP4.5, the rate of mass loss is almost constant over the entire 300 years with a total volume loss equivalent to 20.1 cm sea-20 level increase. A typical thinning pattern for RCP4.5 shows extensive margin thinning and inland retreat of calving fronts after 300 years (Fig. 8). Mass loss near the margin is partially balanced by increased snow accumulation and thickening in the interior.
On the centennial and the three century time scale, the sensitivity of these mass loss projections to the parameter choice in the ice flow model is not very large. For 25 the seven parameter sets that fulfil the calibration criteria best (Table 1), the future sea-level contribution lies within a range of 4 % of the ice loss of the reference model (i.e. ±2 or ±12 mm by 2100 or 2300, respectively). Assuming a 50 % uncertainty in the parameter in the relation for ocean warming-induced discharge increase (Eq. 3), the centennial sea-level contribution deviates by not more than 15 % from the reference model value (with forcing from MIROC-ESM-CHEM). The climate model uncertainty, however, introduces a spread, which is ten times as large for single RCPs (Table 4). Considering the additional scenario spread, the sensitivity to the considered range of ice-sheet model parameters is small.

5
In all climate scenarios, oceanic warming causes additional mass loss from the ice sheet by 2100 AD (upper dark blue columns in Fig. 9). This diagnostics not only comprises the directly induced changes in ice discharge but also their effect on the SMB via geometric adjustments. For individual AOGCM projections, the inclusion of oceanic forcing can explain more than 50 % of the total contribution to sea-level rise by a given 10 time period with an average increase of the total mass loss by ∼ 40 %. In absolute terms, the ocean-induced contribution to sea-level change ranges from 1.8 to 2.6 cm (scenario averages) and 1.1 to 3.2 (full spread) after one century, and from 3.8 to 5.4 cm after three centuries (full spread is 2.3 to 7.4 cm). The oceanic influence on the total ice loss becomes relatively less important for more intense atmospheric warming. 15 It explains about half of the mass loss for RCP2.6 while the share is reduced to 27 % for RCP8.5. This indicates that decreasing SMB and increasing discharge are mutually competitive processes for ice removal at the marine margin. In addition, ice further upstream is efficiently removed by ablation before it actually reaches the marine margin for calving. The oceanic forcing typically induces a diffusive thinning wave at the 20 marine margin which is gradually transmitted inland (Fig. 10a). In areas with a marine margin, this additional thinning wave explains a large share of the total thinning including surface melting under atmospheric warming (Figs. 8 and 10).
In Fig. 9, we also attribute simulated mass changes to either changes in ice discharge, arising from oceanic forcing and inland ice dynamics, or from changes of the 25 mass balance at the ice sheet surface or base (although in all cases, basal melting contributes less than 3 % of the total land ice loss). While increased discharge explains about 40 % of the average mass loss between 2000 and 2010 (light blue columns), its relative contribution generally decreases afterwards and changes in SMB become the dominant factor for mass loss. This is because total ice export across calving fronts eventually falls below year 2000 levels, despite warmer ocean temperatures. Limitations on the ice discharge increase are a direct result of gradual thinning at the marine margins with a fast adjustment of the ice inflow from upstream  but are also a consequence of a retreat of the ice sheet margin back on land. For the 5 CanESM2 model under RCP4.5, the ice sheet loses more than half of its contact area with the ocean by 2300 (Fig. 8). In general, ice discharge increase is more relevant for the total mass loss in scenarios with higher mitigation efforts (RCP2.6, RCP4.5).
The reason is that an ice discharge increase also causes dynamic thinning inland and thereby intensifies surface melting. Surface melting in turn competes with the discharge 10 increase by removing ice before it reaches the marine margin. Margin thinning and retreat limit the ice discharge and increase the relative importance of surface melting in the future volume evolution. The total 2100 ice loss, from SMB changes only, increases by more than 70 % when including ice-ocean interaction. This share is about 42 % of the combined total ice loss in 2100 ( Fig. 9b) but only 10 % of it are directly caused by 15 ice discharge increase at the marine margin. By 2300 AD, the cumulative effect from ice discharge changes are even negative as ice discharge rates have on average fallen below the pre-2000 level between 2000 and 2300.
Detailed flow-line projections of the ice discharge evolution of four major outlet glaciers on Greenland show a general increase by 2100 and 2200 AD (Nick et al.,20 2013). Such a widespread increase of ice discharge is not confirmed by our projections. The glaciers in the Nick et al. (2013) study are however driven with only one specific climate model and only represent the response of four individual, well-studied outlet glaciers. In our large-scale model approach, ice discharge of main outlet glaciers can also show a significant increase while the ice sheet-wide discharge increase is more 25 moderate. In other words, as it occurs many of the smaller outlet glaciers quickly become land-based. Therefore, scaling up the discharge response of only those glaciers with the most prolific ice export is not necessarily representative for the future icedynamic evolution of an entire ice sheet. A generalisation of the discharge evolution of the four outlet glaciers modelled in Nick et al. (2013) to the entire ice sheet is in line with our finding that the relative importance of ice discharge changes to the future ice loss is self-limited by thinning and retreat of ice in contact with the ocean .
In all experiments, the additional effect of basal lubrication on total mass loss is very 5 small, corresponding to an additional sea-level contribution of less than 1 % (Fig. 10b). This is in agreement with results from a parametric approach to link runoff to basal lubrication (Shannon et al., 2013). Lubrication-induced speedup displaces inland ice mass but does in general not remove it. In the upper ablation area, the ice thins as it accelerates, while for melt rates exceeding 2 m yr −1 near the margin, the relative 10 speedup decreases under warming, causing a relative thickening (cf. Eq. 2 and Fig. 1). The reason is that when meltwater export rates exceed a threshold, a channellisation of the basal drainage system is assumed with concurrent reduction of basal lubrication. Ice flow is mainly enhanced close to the equilibrium line where runoff rates cause maximal speedup. This even leads to a negative feedback as the relative thickening 15 of the ablation zone reduces runoff rates through the height-mass balance feedback (Huybrechts et al., 2002).

Summary and conclusion
In this study, we included more dynamic processes in a thermo-mechanically-coupled, three-dimensional ice flow model with the aim to better assess the impact of dynamic 20 changes on the future sea-level contribution of the Greenland ice sheet. For this purpose, climate anomalies were taken from a suite of ten Atmosphere-Ocean General Circulation Models (Table B1)  higher-order ice-dynamic model component allows for a direct inland transmission of the imposed margin perturbations. In this respect, it is in very good agreement with the observational range of 0.7 ± 0.1 mm yr −1 (Shepherd et al., 2012). In this period, increased ice discharge explains ∼ 40 % of the total mass loss of the ice sheet, also in agreement with observational evidence. To our knowledge, this is the first model 5 approach that achieves to reproduce the recent mass loss and explains it by changes in SMB and in ice discharge. The latter changes are attributed to oceanic warming in the surrounding ocean basins. The mean volume loss for the CMIP5 ensemble is however biased low with 0.32 mm yr −1 . This bias arises from the spread in climate models that are not all expected to simulate the observed trend over such a short time pe-10 riod. But the range of mass loss from these CMIP climate models has a maximum at 0.71 mm yr −1 and thereby covers values inferred from observations. For the climate model ensemble, increased ice discharge also explains ∼ 40 % of the total mass loss during the last decade. Accounting for the four RCP scenarios, we find a Greenland ice sheet contribution 15 to global sea-level rise of between 1.4 and 16.6 cm by 2100 AD. For the two lowimpact scenarios, ice loss attains respectively 11.1 and 32.0 cm by 2300 AD. Despite an average increase in mass loss of ∼ 40 % in 2100, when accounting for ice-ocean interaction, mass loss is predominantly caused by changes in SMB. The reason is that ice discharge is limited by margin thinning and retreat but also by a competition with 20 surface melting that removes ice before it reaches the calving front. These geometric limits on ice discharge explains that changes in SMB cause most of the mass loss by 2100. Beyond 2100, modelled ice discharge rates fall below the pre-2000 level and this decrease is compensated by the dominant changes in SMB. The results therefore suggest that the largest source of uncertainty in future mass loss arises from the SMB 25 and the underlying climate change projections, and not from ice dynamics. Our results have implications for attempts to estimate the role of ice discharge on the future mass loss of the Greenland ice sheet. Observed rates of change over the last decade cannot simply be extrapolated over the 21st century on account of a different balance of processes causing mass loss over time. Extrapolating recently observed mass trend changes to a century time scale (Rignot et al., 2011) or linking observed Greenland sea-level trends to temperature change (Rahmstorf, 2007) implies continued glacier acceleration and a multifold increase of the ice discharge (Pfeffer et al., 2008) that is not found attainable in numerical ice-sheet models. Ice discharge at calv-5 ing fronts is self-limited by ice dynamics, supporting the view that centennial mass changes are dominantly driven by SMB changes, and thus by changes in surface climate conditions.

Appendix A: Climate conditions over the last glacial cycle
Temperature history 10 The model spin-up over several glacial cycles requires information on the past climate, which is reconstructed from ice core data. The glacial temperature forcing is obtained from synthesised isotope records representative for central Greenland conditions. For the period prior to 122.6 kyr BP, the forcing reconstruction is based on a synthesised Greenland δ 18 O record derived from Antarctica Dome C using a bipolar seesaw model 15 (Barker et al., 2011). Subsequently, the NGRIP δ 18 O record (Andersen et al., 2004) is used before switching to GRIP information at 103.8 kyr BP (Dansgaard et al., 1993). For the last 4000 years, a direct reconstruction of snow temperatures is available based on a a δ 15 N/δ 40 Ar record from GISP2 (Kobashi et al., 2011  . Thereafter, the temperature reconstruction shows 5 a mismatch of 0.4 • C with the GRIP reconstruction at 4000 years BP. Before splicing these two records, the Kobashi et al. (2011) temperatures are lowered over time with a linear function that removes the past mismatch but keeps the present-day values (Fig. A1). In a final step, the temperature reconstruction is linearly interpolated on time intervals of ten years. 10 Assembling the forcing record in this way prolongs any records exclusively based on Greenland ice cores by several hundred millennia. In addition, the intermediate switch to the NGRIP record gives more reliable information during the late Eemian period than GRIP. This is because of known disturbances in the lower parts of the GRIP ice core prior to 105 kyr BP. The last splice with surface snow temperature reconstructions 15 at GISP2 seems favourable because this reconstruction method was validated against observations and model reconstructions starting in 1850 AD. One remarkable feature of our assembled temperature forcing record is the Little Ice Age cooling on the Greenland ice sheet (Fig. A1). This cold period 200-500 years ago influences our spin-up into the present-day, and causes ice sheet growth up to the beginning of the 20th century.

Appendix B: Breakdown of projections by climate model
For most of the climate model ensemble members (Table B1), air temperature anomalies correlate better with the centennial contribution of the Greenland ice sheet to sealevel change than ocean temperature anomalies. Linear correlation coefficients for air temperature in general exceed 0.7 while this threshold is not surpassed for ocean tem- 25 peratures except in RCP8.5. By 2300, the correlation with ocean forcing dominates for RCP2.6.
3876 The spread in centennial sea-level contributions and atmospheric warming (Fig. B1) reflects both uncertainties in the realised future scenario and differences in the respective AOGCM. Up to 2100 AD, this spread is explained by differences in individual AOGCM projections rather than scenario differences. In particular the three low impact scenarios show a large overlap in AOGCM realisations. By 2300, the spread introduced TCD Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version