The GRISLI-LSCE contribution to ISMIP6, Part 1: projections of the Greenland ice sheet evolution by the end of the 21st century

Polar amplification will result in amplified temperature changes in the Arctic with respect to the rest of the globe making the Greenland ice sheet particularly vulnerable to global warming. While the ice sheet has been showing an increase mass loss in the past decades, its contribution to global sea level rise in the future is of primary importance since it is at present the largest single source contribution behind the thermosteric contribution. The question of the fate of the Greenland and 5 Antarctic ice sheets for the next century has recently gathered various ice sheet models in a common framework within the Ice Sheet Model Intercomparison Project for CMIP6. While in a companion paper we present the GRISLI-LSCE contribution to ISMIP6-Antarctica, we present here the GRISLI-LSCE contribution to ISMIP6-Greenland. We show an important spread in the simulated Greenland ice loss in the future depending on the climate forcing used. The contribution of the ice sheet to global sea level rise in 2100 can be thus as low as 20 mmSLE to as high as 160 mmSLE. The CMIP6 models produce much 10 larger ice sheet retreat than their CMIP5 counterparts. Low emission scenarios in the future drastically reduce the ice mass loss. The mass loss is mostly driven by atmospheric warming and associated ablation at the ice sheet margin while oceanic forcing contributes to about 10 mmSLE in 2100 in our simulations.

fate of the ice sheets in the future. First, in order to save computing time, most of the ice sheet models use various asymptotic approximations (e.g. the shallow ice and shallow shelf approximations or higher-order models) even if more recently models that account explicitly for all the stress components of the Stokes equation at the ice sheet scale have emerged (e.g. Seddik et al., 2012). This difference in term of ice sheet model complexity is a source of uncertainty for future projections. Second, ice sheets respond to a wide spectrum of timescales, from sub-annual to multi-millenial. As a result, diverse methodologies 5 to initialise the models for projection purposes have been developed. For Greenland ice sheet models, these differences in methodologies lead to an even larger uncertainty for future projections than model complexity and explain most of the multimodel spread (Goelzer et al., 2018). A last source of uncertainty lies in poorly known processes, such as sub-glacial processes, or processes that are not included in models due to their complexity or too fine spatial scale, such as outlet glacier dynamics or fracturing. Large international intercomparison exercises are a useful way to quantify these different uncertainties and to infer 10 robust sea level projections into the future.
The Ice Sheet Model Intercomparison Project for CMIP6 (ISMIP6, Nowicki et al., 2016), endorsed by the Coupled Model Intercomparison Project -phase 6 (CMIP6), aims at investigating the role of dynamic Greenland and Antarctic ice sheets in the climate system and to reduce the uncertainty in ice sheet contribution to global sea level rise in the future. Within this 15 framework, stand-alone ice sheet model experiments have recently been carried out by world-wide research groups. Many model experiments using both CMIP5 and CMIP6 climate forcing scenarios until 2100 were conducted with ice sheet models spanning a range of model complexities and using different initialisation techniques. To date, this is the most ambitious intercomparison exercise dedicated to the fate of the Greenland and Antarctic ice sheets in the future. At the Laboratoire des Sciences du Climat et de l'Environnement (LSCE), we participated to this stand-alone intercomparison with the GRISLI model 20 (Quiquet et al., 2018). This model uses the shallow ice and shallow shelf approximations and is relatively inexpensive in term of computational cost. We were thus able to perform all the different experiments of ISMIP6. The aim of this paper is to discuss the role of the forcing uncertainties for future projections of the Greenland ice sheet contribution to global sea level rise when using our model. This individual model response can be put in perspective with respect to the multi-model spread discussed in Goelzer et al. (2020). A companion paper (Quiquet and Dumas, submitted) describes the results for the Antarctic ice sheet. 25 In Sec. 2 we describe briefly the GRISLI ice sheet model as well as the procedure used for its initialisation. We also provide information on the ISMIP6-Greenland forcing methodology and we provide an overview of the different experiments performed.
In Sec. 3 we discuss the results for the different experiments in term of geometry and dynamical changes. We discuss these results in a broader context in Sec. 4 and we conclude in Sec. 5. 30 2 https://doi.org/10.5194/tc-2020-139 Preprint. Discussion started: 2 June 2020 c Author(s) 2020. CC BY 4.0 License.

Model and initialisation
For this work, we use the GRISLI ice sheet model. The model is a 3D thermo-mechanically coupled ice sheet model that solve the mass and momentum conservation equations. The model is fully described in Quiquet et al. (2018) and we only provide here a brief overview of its characteristics.

5
Assuming incompressibility, the mass conservation equation for a grid element is: with H the local ice thickness, BM the mass balance andŪ the vertically averaged horizontal velocity vector. ∇ Ū H is thus the ice flux divergence.
The Stokes momentum equation is solved using asymptotic zero-order approximations. For the whole domain, we assume that 10 the total velocity is simply to superposition of the two main approximations: the shallow ice approximation (SIA) in which the horizontal derivatives of the stress tensor are neglected (vertical shear driven deformation) and the shallow shelf approximation (SSA) in which the vertical derivatives of the stress tensor are neglected (longitudinal stresses predominant). Practically, this means that we use the SSA equation as a sliding law (Bueler and Brown, 2009). Grounded cold base and floating shelves are special cases for which there is infinite, respectively none, friction at the base. Elsewhere, friction is assumed to follow a 15 Weertman (1957) power law with a till that allows viscous deformation: where τ b is the basal drag, β is the basal drag coefficient and U b is the basal velocity. The basal drag coefficient is spatially variable but constant in time (except in specific cases such as during the inversion procedure). 20 Similarly to what has been done with GRISLI for the initMIP-Greenland experiments (Goelzer et al., 2018), we used here an inverse procedure to initialise the model at the start of the historical experiment. We mostly followed the iterative method of Le clec'h et al. (2019b) which consists at yielding the map of the basal drag coefficient β that minimises the ice thickness error with respect to observations. To this aim, we first run a 30 kyr experiment with fixed topography and perpetual present-day climate forcing in order to compute the thermal state of the ice sheet in agreement with the boundary conditions. From this, we 25 do multiple 200-yr long experiments under constant present-day climate forcing but with an evolving topography. During the first 20 years of these experiments, we adjust the basal drag coefficient to minimise the ice thickness mismatch with respect to the observations. Each iteration starts from the exact same initial condition, except that the basal drag coefficient is different.
Also, the ice thickness error at the end of the 200-yr long experiment is used to facilitate convergence towards the observed ice thickness through a local basal drag modification. This modification on the basal drag consist in finding an ice flux on the 30 simulated topography as close as possible to the balance ice flux on the observed topography. After a few 200-yr experiments, we repeat the thermal equilibrium computation so that it accounts for the newly inferred basal drag coefficient. For this work we performed more than ten thermal equilibrium experiments, each one followed by five iterations of 200 years.
At the end of the iterative process, we use the last inferred basal drag coefficient together with the corresponding thermal state to run a short relaxation experiment of 20 years. The end of this relaxation experiment defines our initial state which is used to begin the historical experiment hist and the control experiment ctrl (see Sec. 2.3).

5
Our ice thickness and bedrock topography of reference is the BedMachine v.3 (Morlighem et al., 2017). This dataset is used as a target for our iterative procedure to infer the basal drag coefficient. It is also used as the starting topography for the short relaxation that defines our initial state. Our present-day climate forcing, namely annual near-surface air temperature and annual surface mass balance, comes from the MAR v3.9 (Fettweis et al., 2013(Fettweis et al., , 2017 forced at its boundary by MIROC5, averaged 10 over the 1994-2015 period. On top of this climate forcing, we also add a strongly negative surface mass balance term of -15 m yr -1 outside the present-day ice mask in the observational dataset in order to avoid inconsistencies between the climate forcing and the initial ice sheet geometry. The model is run on a Cartesian grid at 5 km resolution covering the Greenland ice sheet using a stereographic projection.

15
The ISMIP6-Greenland working group distributed atmospheric and oceanic forcings to drive individual ice sheet models. They also suggest a forcing methodology so that participating models are run using a common framework. Full description of the methodology is available in Nowicki et al. (2020) and only a summary is presented here.
For the atmospheric forcing, MAR v3.9 has been run from 1950 until 2100 forced at its boundaries by a selection of CMIP5 20 and CMIP6 general circulation model (GCM) outputs. To force the ice sheet models, yearly anomalies of near-surface air temperature and surface mass balance are provided. These anomalies were constructed as the difference of a given yearly value with the climatology over the reference period 1960-1989. In addition, to account for the surface elevation feedback on temperature and surface mass balance, yearly values of atmospheric gradients for these two variables are also provided. These  Ice-ocean interactions for the Greenland ice sheet is most of the time poorly represented amongst ISMIP6-Greenland participating models. This is mostly due to the fact that the spatial scale needed to represent such interactions is out of reach for most models. This is also the case for GRISLI, where the 5 km resolution grid is too coarse to capture marine-terminating outlet glaciers. To cope with this problem, retreat masks for outlet glaciers have been available in ISMIP6-Greenland. They were 30 obtained with simple parametrisations calibrated and tested against observational datasets (Slater et al., 2019). These masks provide, for a given resolution, the fraction of the grid that becomes ice free and they are used to impose a specific retreat rate of the marine front. For each climate forcing, three retreat masks are available for different oceanic sensitivities (low, medium and high). Since our model does not account for partially glaciated grid cell, the fractional information given by the retreat masks is used to reduce the local ice thickness with respect to a reference ice thickness (i.e. the ice thickness evolution for the outlet glaciers is imposed). The reference ice thickness could have been chosen as the ice thickness at a specific time (e.g. the ice thickness at the end of the historical experiment). However, in doing so, we may create strong discontinuities in ice thickness when the retreat mask is used for the first time. For this reason, we choose instead the value of the local ice thickness at the time when the imposed retreat starts to play as a reference ice thickness. to CMIP6, sensitivity to the greenhouse gas emissions scenarios, the respective role of oceanic forcing with respect to atmospheric forcing, and to quantify the uncertainty regarding the oceanic forcing.
The core experiments consist (Tier 1) in a selection of three CMIP5 climate models (MIROC5, NorESM and HadGEM2-ES) run under the RCP8.5 scenario for greenhouse gases. In addition, MIROC5 was chosen to be run with a different RCP scenario (RCP2.6) and using different oceanic sensitivities (high and low in addition to medium). Tier 2 has two subsets: an extended 20 ensemble with three additional CMIP5 models using RCP8.5 and an other with four CMIP6 models. Amongst CMIP6 models, CNRM-CM6 has been run under two scenarios: a pessimistic (SSP585) and an optimistic (SSP126) scenario. Tier 3 has also two subsets. The first one aims at quantifying the respective role of the oceanic forcing with respect to atmospheric forcing, running the ice sheet models only with one of this forcing at a time. Three climate models were selected (MIROC5, CSIRO-Mk3.6 and NorESM) and as in Tier 1, MIROC5 was run for two greenhouse gases scenarios and different oceanic sensitivities.

25
Finally, the second subset of Tier 3 contains the ten climate models (CMIP5 and CMIP6) each time run with the two additional oceanic sensitivities (high and low). CNRM-CM6 is the only one in this subset that has run under two emission scenario (SSP585 and SSP126). We performed all these experiments with the GRISLI ice sheet model.
In addition to the these projection experiments, we also perform two control experiments in which the climate forcing remains 30 unchanged, being our reference climate forcing used during the initialisation procedure (zero anomaly). The control experiment ctrl starts from the initial state resulting from our initialisation procedure and covers the 1995-2100 period (106 years). The ctrl_proj experiment starts in January 2015, alike the projection experiments, and runs for 86 years under a constant climate forcing.
We aim here at providing a detailed description of the historical experiment hist and the model response under the various forcings of the projection experiments. While some information is given in this section, the reader is invited to refer to Goelzer et al. (2020) to compare in details the response of GRISLI to other participating models.

Present-day simulated ice sheet 5
At the end of the historical experiments hist, with a value smaller than 30 m, GRISLI shows the lowest ice thickness root mean squared error (RMSE) with respect to the observations of Morlighem et al. (2017) amongst the ISMIP6-Greenland participating models (Goelzer et al., 2020). This is a result of the initialisation procedure we use (Sec. 2.1) that includes only a short relaxation of 20 years. With an historical experiment of 20 years only, the model has no time to depart strongly from the observations. The map of the ice thickness difference with respect to observations is shown in Fig. 1a. The model shows a very good 10 agreement with the observations for most of the ice sheet, except at specific locations at the margin. In particular, South-East Greenland is the least well reproduced with local errors greater than 200 metres. In the region of Kangerdlugssuaq and Helheim glaciers, there is an ice thickness overestimation near the glacier termini and an underestimation upstream. These differences with the observations can be due to the fact that this area is particularly difficult to model since it has a complex surface mass balance pattern with very strong horizontal gradients and also a rough topography that is not necessary well captured at 5 km 15 resolution.
Some of the ice thickness mismatch with respect to the observations can be partly explained by error related to ice dynamics.
GRISLI has indeed an ice velocity RMSE with respect to the observations (Joughin et al., 2016) of about 35 m yr -1 , making the model the sixth worst model out of 21 (Goelzer et al., 2020). Our initialisation procedure favours a good match of the 20 simulated ice thickness with respect to observations but it does not include any constraints on the ice velocity. It is thus not particularly surprising that GRISLI performs best in term of ice thickness than in term of ice velocity. Since ice velocity is a very heterogeneous variable, it is sometimes convenient to use the logarithm of the velocity instead of the absolute velocity.
In doing so, the RMSE is about 0.55 (eleventh worst value out of 21). In using the logarithm of the velocity, GRISLI slightly improves compared to the other participating models, meaning that the errors are mostly localised in areas of high velocities.
25 Fig. 2a shows the absolute simulated velocity, to be compared to the observations in Fig. 2b. The pattern is generally well reproduced and the model is able to reproduce the localisation of the major existing ice streams. However, the velocity of the ice streams is not always in agreement with the data. The Northern and Western ice streams are generally too slow with an underestimation reaching more than 500 m yr -1 for the Jakobshavn, Petermann and North East Greenland ice stream (NEGIS) glacier termini (Fig. 2c). On the contrary, the South East glaciers, Kangerdlugssuaq and Helheim, are too slow in the model.

30
For the northern and western regions, the errors in ice thickness are small meaning that the ice velocity mismatch there cannot be reduced within our initialisation procedure which only minimises the ice thickness error. This is somewhat different for the South Eastern region, where there are important errors in ice thickness. However, there is an important positive bias in ice thickness at the ice sheet margin that tends to produce very high ice flow (very low basal drag coefficient to reduce this bias). Since the SSA equation is elliptic, the low basal drag at the margin has a regional impact on ice flow, which tends to produce an underestimation of the ice thickness further inland. While we strongly overestimate the velocity in this area, the ice thickness at the margin is still overestimated. This suggests that the surface mass balance used in our reference climate is probably overestimated in this region.

5
The ice sheet model drift can be assessed by examining Fig. 1c. The ice thickness drift in the control experiment ctrl_proj is generally very small (lower than 10 metres) with only a few regions with higher values. Here again, the Kangerdlugssuaq and Helheim glacier regions show the largest model drift with a local increase in ice thickness of more than 100 metres near the glacier termini. Overall the ice volume drift is negligible over the duration of the control experiment (86 years), also 10 because of some compensating biases (see also Fig. 3). In addition to a the ice thickness drift, the model simulates a drift in velocities during the duration of the control experiment (Fig. 2d). For most of the ice sheet the velocity change is small and only reaches more than 1 m yr -1 at the ice sheet margins. The largest changes concern the glaciers in South-East Greenland such as the Helheim and the Kangerdlugssuaq glaciers where locally, at the termini, there can be an acceleration by more than 1000 m yr -1 .

25
CMIP6 models show generally a much larger climate sensitivity than their equivalent in the former CMIP5 generation (Forster et al., 2020). This has important consequences on the projected Greenland ice sheet. The ice volume evolution for the four CMIP6 models under SSP585 scenario is shown in Fig. 4. The CMIP6 models produce systematically higher ice loss than the CMIP5 models. The two most sensitive CMIP6 models (UKESM1-CM6 and CESM2) almost double the ice loss with respect 30 to the most sensitive CMIP5 models (IPSL-CM5-MR and MIROC5). The ice loss thus reaches -60000 km 3 (140 mmSLE) by the end of the century.
Two climate models have been run under two scenarios for the evolution of future atmospheric greenhouse gases. The ice loss for the two scenarios of the climate models is shown in Fig. 5. The CMIP5 (MIROC5) and CMIP6 (CNRM-CM6) responses to the change in greenhouse gas scenario (RCP8.5 to RCP2.6 and SSP585 to SSP126 respectively) is very similar. There are very small differences for the first half of the century but after 2060 the pessimistic scenario produces substantial additional mass loss with respect to the optimistic scenarios. By the end of the century, the pessimistic greenhouse gas scenario produces 5 roughly -25000 km 3 (55 mmSLE) of additional ice loss with respect to the optimistic scenario. The future evolution of atmospheric greenhouse gas mixing ratio is thus a major driver for the Greenland ice mass at the century time scale.
The spatial pattern of ice loss by the end of this century is shown in Fig. 6. For this figure we have chosen four projection experiments that show contrasted integrated ice mass loss by 2100: the CISRO-Mk3.6 under RCP8.5 which produces a small 10 integrated ice loss (Fig. 6a), the MIROC5 under RCP8.5 which produces an important mass loss (Fig. 6b), the MIROC5 under RCP2.6 to show the impact of low emission scenario (Fig. 6c) and UKESM-CM6 under SSP585 with a high oceanic sensitivity which produces the highest mass loss amongst all the different experiments (Fig. 6d). While the amplitude of ice thickness change is drastically different amongst these experiments, the spatial pattern is similar. The major signal is a substantial widespread ice thickness decrease at the margin of the ice sheet. If the ice thickness decrease is about 50 m for the least

Importance of the oceanic forcing
The uncertainty that arises from the oceanic forcing can be evaluated thanks to the different glacier retreat scenarios (low, medium and high sensitivity to oceanic forcing). In Fig. 3 is represented on the right-hand side the uncertainty that arises from 25 the oceanic forcing for the individual CMIP5 models. In 2100 the ice volume loss difference between the low and high oceanic sensitivities is generally of about -5000 km 3 (less than 10 mmSLE). Without being negligible, the oceanic sensitivity for a given climate scenario is nonetheless relatively small compared to the spread amongst the different CMIP5 climate models used. For the CMIP6 experiments, the uncertainty that comes from the oceanic forcing is almost doubled with respect to the CMIP5 experiments, with about 10000 km 3 (~20 mmSLE) of ice loss difference from the low to high oceanic sensitivity (Fig. 4) but 30 these CMIP6 models also produce much greater ice loss.
We also performed experiments in which we isolate the response of the model that arises from the atmospheric forcing only only the atmospheric perturbation is taken into account. Conversely, for the ocean only (OO) experiments, there is no atmospheric perturbation (as in the control ctrl_hist experiment) but we do impose a retreat rate for the outlet glaciers. The ice volume evolution for these experiments is shown in Fig. 7. The OO experiments produce almost identical ice volume evolutions amongst the different GCMs. This means that even if the glacier retreat is subject to uncertainties, with the methodology of Slater et al. (2019) it is nonetheless only weakly sensitive to the differences in the climate forcing used to elaborate it.
5 Fig. 7 also shows that the atmospheric forcing is the main driver for ice loss for the GCMs that produce important ice loss.
Also, the sum of the ice loss of AO and OO experiments approximate closely the ice loss simulated when using the full forcing.

Change in ice dynamics
Climate forcing, and its associated ice sheet geometry change, leads to a change in the dynamics of the Greenland ice sheet.  Fig. 8b shows the difference in ice flux convergence in 2100 for a selected climate forcing with respect to the control ctrl_proj experiment. This can be considered as the dynamical contribution to ice thickness change. The pattern mostly 20 follows the one of velocity change (Fig. 8a). There is an important ice flux convergence at the margins that tends to partially compensate the decrease in surface mass balance. Conversely, upstream regions show a slightly negative convergence. This pattern is similar amongst the different climate forcings.

25
In order to minimise the initial error in ice thickness with respect to the observations, we have used an inverse procedure that optimally tune the basal drag coefficient. In doing so, we produce a simulated ice sheet that is in quasi-equilibrium with the climate forcing (minimal ice thickness drift). In reality, the Greenland ice sheet is far from being at equilibrium with present-day climate since it has been loosing ice at an accelerated rate over the last four decades .
This means that, by constructions, our simulations underestimate the Greenland ice sheet contribution to future sea level rise.

30
Some promising alternatives exist, for example using data assimilation of observed velocities in a transient ice sheet simulation (Gillet-Chaulet, 2020). These methods require however a complex data assimilation framework, currently not implemented in our ice sheet model. Instead, we plan to modify the inverse procedure of Le clec'h et al. (2019b) by incorporating the ice thickness change inferred by gravimetry as an additional constraint in order to improve on the initial state of the Greenland ice sheet.
In addition, the inferred basal drag coefficient during the initialisation procedure is left unchanged for the duration of the historical and projection experiments. This is probably an important and unjustified approximation since the basal conditions are 5 susceptible to respond to changes in ice geometry and, eventually, basal hydrology. To assess the importance of basal drag coefficient changes for our projections, we perform a new set of experiments using the MIROC5 climate forcing under RCP8.5 with a medium oceanic sensitivity. For these simulations, we apply a spatially uniform modification factor to reduce or increase the value of the basal drag coefficient after the year 2045. The ice volume difference in 2100 with respect to the experiment with no modification of the basal drag coefficient is shown in Fig. 9a,b. The ice volume change in response to small perturbations 10 of the basal drag coefficient is relatively linear and limited. Thus, a perturbation of 20% results in less than 5000 km 3 volume change, which translates to less than 10 mmSLE. This means that it is unlikely that basal condition changes in the future could produce a drastically different ice volume change in 2100. This also suggests that a slightly different basal drag coefficient inferred during our initialisation procedure will produce a similar ice volume evolution in the projection experiments. We repeat this kind of sensitivity experiment for the enhancement factor (Fig. 9c,d). Here again, the response in term of ice volume is 15 small and relatively linear.
One additional limitation of the inverse procedure is that it does not take into account the impact of the last glacial cycle on ice temperatures. Our internal temperature field is the result of a long thermo-mechanical equilibrium under perpetual present-day forcing and as such, it is necessarily overestimated. In addition to an underestimated ice viscosity, this has also consequences on 20 the simulated basal temperature and, as a result, on the regions where sliding occurs. This might affect the dynamical response of the model to future climate change. If earlier studies have already identified these limitations (e.g. Rogozhina et al., 2011;Yan et al., 2013;Seroussi et al., 2013;Le clec'h et al., 2019b), our inverse procedure does not allow for a quantification of these limitations. Given the relatively low computational cost of GRISLI, one alternative would be to perform muti-millenial palaeo integrations to infer the initial state used for the projections. This strategy generally leads to a larger ice thickness error 25 with respect to present-day observations but has the advantage to have a thermal state consistent with the model physics and with the palaeo temperatures. While the ISMIP6-Greenland participating models either choose one or the other initialisation technique (Goelzer et al., 2020), it would be very informative to have two drastically different initialisation methods for a given ice sheet model.

30
Finally, the forcing methodology used for ISMIP6-Greenland does account for the vertical elevation feedback on temperature and surface mass balance. However, other feedbacks are at play for the future evolution of the Greenland ice sheet. Notably, the MAR model used to compute the forcing fields does not account for topography and ice mask changes. The effect of these changes is probably limited for moderate ice sheet retreat (Le clec'h et al., 2019a). However, since CMIP6 models produce a much greater retreat, they could also induce more important feedbacks if MAR was bi-directionally coupled to an ice sheet 35 model. In addition, the effect of Greenland ice loss on the ocean is also not taken into account with the forcing methodology followed. While, the oceanic forcing seems not to be the major driver for future Greenland ice loss, glacier retreat in the future should ideally take into account the oceanic circulation changes in the fjords related to freshwater discharge from ice sheet melting.

Conclusions
In this paper we have presented the GRISLI-LSCE contribution to ISMIP6-Greenland. Independently from the climate forcing used to drive the ice sheet model, we have shown that the Greenland ice sheet is systematically loosing ice in the future.
However, the magnitude of the mass loss by 2100 is very sensitive to the climate forcing. Under a business as usual scenario for the greenhouse gas emission (RCP8.5 or SSP585), the mass loss translates into a Greenland ice sheet contribution to global 10 sea level rise that ranges from 35 to 160 mmSLE. However, with an optimistic greenhouse gas emission scenario (RCP2.6 or SSP126) the mass loss can be significantly reduced. The CMIP6 models tend to produce much larger ice loss due to their higher climate sensitivity with respect to the one of the CMIP5 models. While the oceanic forcing contributes to ice loss by about 10 mmSLE in 2100, the Greenland ice mass loss in the future is mostly driven by a larger ablation at the ice sheet margin in the future. This suggests that this process should be carefully implemented in ice sheet models aiming at simulating the 15 Greenland ice sheet evolution at the century scale.

Data availability
The GRISLI outputs from the experiments described in this paper are available on the zenodo repository with digital object identifier 10.5281/zenodo.3784665. The outputs in the zenodo repository are the standard GRISLI outputs on the native 5 km grid and, as a result, they may slightly differ from the post-processed outputs available on the official CMIP6 archive on the 20 Earth System Grid Federation (ESGF). In order to document CMIP6's scientific impact and enable ongoing support of CMIP, users are obligated to acknowledge CMIP6, the participating modelling groups, and the ESGF centres (see details on the CMIP    and ice volume contributing to sea level rise ((b) and (d)).