Articles | Volume 14, issue 10
Research article
02 Oct 2020
Research article |  | 02 Oct 2020

Sensitivity of Greenland ice sheet projections to spatial resolution in higher-order simulations: the Alfred Wegener Institute (AWI) contribution to ISMIP6 Greenland using the Ice-sheet and Sea-level System Model (ISSM)

Martin Rückamp, Heiko Goelzer, and Angelika Humbert

Projections of the contribution of the Greenland ice sheet to future sea-level rise include uncertainties primarily due to the imposed climate forcing and the initial state of the ice sheet model. Several state-of-the-art ice flow models are currently being employed on various grid resolutions to estimate future mass changes in the framework of the Ice Sheet Model Intercomparison Project for CMIP6 (ISMIP6). Here we investigate the sensitivity to grid resolution of centennial sea-level contributions from the Greenland ice sheet and study the mechanism at play. We employ the finite-element higher-order Ice-sheet and Sea-level System Model (ISSM) and conduct experiments with four different horizontal resolutions, namely 4, 2, 1 and 0.75 km. We run the simulation based on the ISMIP6 core climate forcing from the MIROC5 global circulation model (GCM) under the high-emission Representative Concentration Pathway (RCP) 8.5 scenario and consider both atmospheric and oceanic forcing in full and separate scenarios. Under the full scenarios, finer simulations unveil up to approximately 5 % more sea-level rise compared to the coarser resolution. The sensitivity depends on the magnitude of outlet glacier retreat, which is implemented as a series of retreat masks following the ISMIP6 protocol. Without imposed retreat under atmosphere-only forcing, the resolution dependency exhibits an opposite behaviour with approximately 5 % more sea-level contribution in the coarser resolution. The sea-level contribution indicates a converging behaviour below a 1 km horizontal resolution. A driving mechanism for differences is the ability to resolve the bedrock topography, which highly controls ice discharge to the ocean. Additionally, thinning and acceleration emerge to propagate further inland in high resolution for many glaciers. A major response mechanism is sliding, with an enhanced feedback on the effective normal pressure at higher resolution leading to a larger increase in sliding speeds under scenarios with outlet glacier retreat.

1 Introduction

Climate change is the major driver of global sea-level rise (SLR), which has been shown to be accelerating (Nerem et al.2018; Shepherd et al.2019). The Greenland ice sheet (GrIS) has contributed to about 20 % of sea-level rise during the last decade (Rietbroek et al.2016). Holding in total an ice mass of  7.42 m sea-level equivalent (SLE; Morlighem et al.2017), its future contribution poses a major societal challenge. Since 1992, the GrIS mass loss has been controlled on average 52 % by surface mass balance (SMB), with the remainder of 48 % being due to increased ice discharge of outlet glaciers into the surrounding ocean (Shepherd et al.2019).

While the relative importance of outlet glacier discharge for total GrIS mass loss has decreased since 2001 (Enderlin et al.2014; Mouginot et al.2019) and is expected to decrease further in the future (e.g. Aschwanden et al.2019), it remains an important aspect for projecting future sea-level contributions from the ice sheet on the centennial timescale (Goelzer et al.2013; Fürst et al.2015). A (non-linear) dynamic response of the ice sheet is caused by changes in the atmospheric and oceanic forcing that may trigger glacier acceleration and thinning of outlet glaciers. Moreover, processes such as SMB and ice discharge are mutually competitive in removing mass from the ice sheet (Goelzer et al.2013; Fürst et al.2015). Beside this interplay, a simple extrapolation of observed GrIS mass loss trends over the next century is not justified, as high temporal variations in SMB and glacier acceleration are apparent (e.g. Moon et al.2012). Therefore, reliable ice sheet models (ISMs) forced with future-climate data must be driven for policy-relevant sea-level projections on century timescales.

The Ice Sheet Model Intercomparison Project (ISMIP6; Nowicki et al.2016, 2020a) is an international community effort striving to improve sea-level projections from the Greenland and Antarctic ice sheets. Based on previous efforts like SeaRISE (Bindschadler et al.2013; Nowicki et al.2013) and ice2sea (e.g. Gillet-Chaulet et al.2012), ISMIP6 continues to fully explore the sea-level rise contribution and associated uncertainties. The effort is aligned with the Coupled Model Intercomparison Project Phase 6 (CMIP6; Eyring et al.2016) to provide input for the upcoming assessment report of the Intergovernmental Panel on Climate Change (IPCC AR6). The general strategy is to use outputs from CMIP5 and CMIP6 climate models to derive atmosphere and ocean fields for forcing ISMs. Goelzer et al. (2020a) and Nowicki et al. (2020b) analyse the future sea-level contribution from multi-model ensembles for ISMIP6-Projections-Greenland. The major aim of Goelzer et al. (2020a) is to provide future sea-level change projections and related uncertainty in a consistent framework.

Despite substantial progress in ice sheet modelling in the last few decades and years, a challenging goal remains to narrow uncertainties and improve the reliability of future sea-level projections from the two big ice sheets. To date, it is recognized that the largest uncertainty sources are related to the initialization of the ISM or stem from external forcing (Goelzer et al.2018, 2020a). Goelzer et al. (2018) compared the initialization techniques used by different ice-sheet-modelling groups. The schematic forward experiment was not designed to estimate realistic sea-level contribution, but it provides valuable insights into how the initial state of an ISM affects the ice sheet response. Under a predefined SMB anomaly, mass losses reveal a large spread. Although the spread is attributed to the broad diversity in model approaches, initMIP-Greenland shows notable improvements (e.g. a reduced model drift) and more consistent results compared to earlier large-scale intercomparison exercises.

Interestingly, the estimated sea-level contributions show a dependence on grid resolution (Fig. 1). ISM versions with multiple grid resolutions demonstrate that coarser grid resolutions tend to produce a slightly larger mass loss. This effect is partly due to the methodological approach by considering an SMB anomaly that is based on the present-day observed SMB. That means ISMs with initial areas larger than those observed are subject to more and stronger melting and sharper transitions in SMB. Therefore, coarse-resolution models not rendering the present-day ice margin perfectly will likely overestimate ablation.

Figure 1Results from the initMIP-Greenland exercise (Goelzer et al.2018). Sea-level contribution vs. (minimum) horizontal grid resolution of each participating ISM. Equal model versions but different grid resolutions are connected with a coloured line. ISSM model versions are connected with a coloured dashed line. Note the logarithmic scale of the x axis. For unstructured meshes the finest resolution is displayed.


However, increasing the spatial resolution comes with the ability to resolve the geometry and to track outlet glacier behaviour (Greve and Herzfeld2013; Aschwanden et al.2016). Some previous works have focused on the dependence of future mass loss of the GrIS on grid resolution (Greve and Herzfeld2013; Aschwanden et al.2019). In these studies no clear conclusion on how the resolution affects the mass loss was found. This was partly explained by the competing tendencies of SMB and ice discharge that are differently resolved by the adopted resolutions. A separation of both responses in future-projection experiments would elucidate how these two main sea-level contributors from the GrIS are affected by the horizontal resolution. Most likely coarse grids underestimate ice discharge as ice flow patterns and cross sections of outlet glacier geometries are not well captured (Greve and Herzfeld2013; Aschwanden et al.2016). High-resolution models, in turn, require a larger amount of computational resource. Unfortunately, when increasing the resolution, simple approximations to the momentum balance do not provide an accurate solution (Pattyn et al.2008). This limitation takes place particularly at the ice sheet margin and at outlet glaciers where all terms in the force balance become equally important (e.g. Pattyn and Durand2013). Due to the intensive computational resources needed to solve the full-Stokes equation, higher-order (HO) approximations provide a good compromise to balance model accuracy and computational costs on centennial timescales.

Determining whether increased model resolution is worth the extra computation time would be valuable to make progress in narrowing uncertainties in ice sheet projections, even if only by a few per cent. The ISMIP6-Projections-Greenland shows that models with low and high resolution are found at the upper and lower bound of sea-level contribution, though no specific analysis of the grid resolution has been performed.

The main intention of this paper is to complement the study by Goelzer et al. (2020a) by evaluating the sensitivity of the simulated GrIS response to global warming due to different horizontal grid resolutions by one single ISM. Besides running the full scenarios (i.e. both oceanic and atmospheric forcing considered), we aim to explore the grid-resolution dependence on atmospheric and ocean forcing separately. Therefore, the full scenarios are complemented with experiments where either a changing SMB or the interaction of the glacier with the ocean is omitted. The simulations of the GrIS are performed with the Ice-sheet and Sea-level System Model (ISSM; Larour et al.2012) and adopted spatial resolutions ranging from medium to high (4 and 0.75 km at fast-flowing outlet glaciers, respectively). The future scenarios build on climate-forcing data from the CMIP5 global circulation model (GCM) MIROC5 under the Representative Concentration Pathway (RCP) 8.5 (Moss et al.2010) following the ISMIP6 protocol (Nowicki et al.2020a).

2 Methods and experiments

Before presenting the concept of this study, we aim to address the terminology used for clarity. Following the glossary in Cogley et al. (2011), ice discharge is computed as the product of ice thickness h and the depth-averaged velocity v¯. In the following, the lower-case q (Gt a−1 m2) refers to the local ice discharge at a point and the upper-case Q (Gt a−1) refers to the glacier-wide quantity (analogous to other quantities such as glacier-wide calving D and local calving d). Quantities at the margin are reckoned to be in the normal direction.

2.1 Ice flow model ISSM

The model applied here is the Ice-sheet and Sea-level System Model (ISSM; Larour et al.2012). It has been applied successfully to the GrIS in the past (Seroussi et al.2013; Goelzer et al.2018; Rückamp et al.2018, 2019a) and is also used for studies of individual drainage basins of Greenland, e.g. the Northeast Greenland Ice Stream (Choi et al.2017), Jakobshavn Isbræ (Bondzio et al.2017) and Petermann Glacier (Rückamp et al.2019b). ISSM is designed to use variable ice flow approximations ranging from shallow-ice approximation to full Stokes and also has the capability to perform inverse modelling to constrain unknown parameters.

Here, we make use of the Blatter–Pattyn approximation (Blatter1995; Pattyn2003) to obtain the most accurate solution even though computational time is increased compared to simpler models (e.g. Aschwanden et al.2019). The system of equations is solved numerically with the finite-element method, and state variables are computed on each vertex of the mesh using piecewise linear finite elements. The ice rheology is treated with a regularized Glen flow law (Glen1955), a temperature-dependent rate factor for cold ice and a water-content-dependent rate factor for temperate ice (Lliboutry and Duval1985).

At the ice base, sliding is allowed everywhere and the basal drag τb follows a linear viscous law (Weertman1957; Budd et al.1984):

(1) τ b , i = - k 2 N v b , i ,

where vb,i is the basal velocity vector in the horizontal plane and i=x,y. Although this type of friction law is often used in ice sheet modelling (Morlighem et al.2010; Price et al.2011; Gillet-Chaulet et al.2012; Seroussi et al.2013), it implies that the basal drag can increase without a bound. It was shown that inducing an upper ratio of τbN (Iken's bound) is more justified (Iken1981; Schoof2005; Gagliardini et al.2007; Leguy et al.2014; Joughin et al.2019). However, we choose this type of friction law as it is commonly used in ISMs making use of inverse methods to constrain the basal friction (Morlighem et al.2010; Larour et al.2012; Seroussi et al.2013; Perego et al.2014; Gladstone et al.2014).

The friction coefficient k2 is assumed to cover bed properties such as bed roughness. The effective pressure is defined as N=ϱigh+min(0,ϱwgzb), where h is the ice thickness; zb is the glacier base; and ϱi=910kgm-3 and ϱw=1028kgm-3 are the densities for ice and seawater, respectively. The parameterization accounts for full water-pressure support from the ocean wherever the ice sheet base is below sea level, even far into its interior where such a drainage system may not exist. At marine-terminating glaciers, water pressure is applied, and zero pressure is applied along land-terminating ice cliffs. A traction-free boundary condition is imposed at the ice–air interface.

The ISSM model domain for the Greenland ice sheet covers the present-day main ice sheet extent and includes the current floating ice tongues (e.g. Petermann, Ryder and 79 N glaciers). The geometric input is BedMachine v3 (Morlighem et al.2017). The thickness, bedrock and ice sheet mask is clipped to exclude glaciers and ice caps surrounding the main ice sheet. The initial ice sheet mask is manually retrieved from the data coverage of the MEaSUREs velocity data set (Joughin et al.2016, 2018) to ensure an available target for the basal-friction inversion employed (Sect. 2.3). A minimum ice thickness of 1 m is applied. Grounding-line evolution is treated with a subgrid parameterization scheme, which tracks the grounding-line position within the element (Seroussi et al.2014). A subgrid parameterization on partially floating elements for basal melt is applied (Seroussi and Morlighem2018). The basal-melt rate below floating tongues is parameterized with a Beckmann–Goosse relationship (Beckmann and Goosse2003). In this parameterization ocean temperature and salinity are set to 1.7 C and 35 Psu, respectively. The melt factor is roughly adjusted such that melting rates correspond to literature values (e.g. Wilson et al.2017; Rückamp et al.2019b). However, at most locations the grounding line coincides with the calving front. Except for the floating-tongue glaciers Petermann, Ryder and 79 N, the subgrid schemes at the grounding line are not applied. The treatment of the calving-front evolution depends on the experimental setup and is explained in Sect. 2.4 and 2.5.2.

Table 1Summary of models and their mesh characteristics. Computational time is based on a projection run under MIROC5 RCP 8.5 and medium ocean forcing.

a Intel Xeon Broadwell CPU E5-2697 v4, 2.3 GHz, on the AWI cluster system Cray CS400.
b Intel Xeon Broadwell CPU E5-2695 v4, 2.1 GHz, on the German Climate Computing Center (DKRZ) cluster system Mistral.

Download Print Version | Download XLSX

As we expect the grounding line not to retreat too excessively in the projections, model calculations with ISSM are performed on a horizontally unstructured grid which remains fixed in time. To limit the number of elements while maximizing the horizontal resolution in regions where physics demands higher accuracy, the horizontal mesh is generated with a higher resolution of REShigh in fast-flowing regions (observed ice velocity > 200 m a−1) and a coarser resolution of RESlow in the interior. This adaptive strategy allows a variable resolution in key areas of the ice sheet, e.g. marine-terminating outlet glaciers. Experiments are carried out at four different horizontal grid resolutions with REShigh equal to 4, 2, 1 and 0.75 km (Table 1 and Fig. 2). The distribution of mesh vertices at numerous outlet glaciers is depicted in Figs. S6 to S19 in the Supplement. In Fig. 3, the interpolated bed elevation for two selected regions and grid resolutions is illustrated. Overall, the bedrock topography of the finer resolution shows deeper and fjord-like troughs, which is closer to the BedMachine data set.

Figure 2Horizontal mesh resolution (km) used for G4000 (a) and G750 (b). Data are clipped at 0.5 and 10 km. The horizontal resolution of a triangle is defined by its minimum edge length. The grey line delineates the initial ice domain. Grey grid lines indicate 100 km. The black boxes indicate the northwestern and southeast subsets used in following figures.

Figure 3G4000 and G750 bedrock topography for the northwestern region (a, b, respectively), and for the southeast region (c, d, respectively). Region subsets are shown in Fig. 2.

Independent of the horizontal resolution, the vertical discretization comprises 15 terrain-following layers, refined towards the base where vertical shearing becomes more important. Please note that G1000 and G750 correspond to the ISMIP6 contributions AWI-ISSM2 and AWI-ISSM3, respectively, in Goelzer et al. (2020a).

2.2 Overview of experiments

The ISMIP6 protocol requests the initialization mode prior to or on the ISMIP6 projection start date. If the initialization date is before the start of the projections, a short historical run is needed to advance the ISM from the reference date to the end of 2014. From this date, the future-climate scenarios branch off. Unforced constant-climate control experiments are defined to capture the model drift with respect to the ISM reference climate and the ISMIP6 projection start date. The set of experiments are described in the following sections but can be summarized as follows:

  • Initialization. Experiment to retrieve the initial state of the model is performed.

  • ctrl. Experiment where the climate is held constant to the reference climate (from January 1991 to end of December 2100) is performed.

  • ctrl_proj. Experiment where the climate is held constant to the reference climate (from January 2015 to end of December 2100) is performed.

  • Historical. Experiment to bring model from the initialization state to the ISMIP6 projection start date (from January 1991 to end of December 2014) is performed.

  • Projection. Future-climate scenario (from January 2015 to end of December 2100) is projected.

2.3 Initialization experiment

The initialization state of ISSM is based on data assimilation and inversion for determining the basal-friction coefficient. Before the inversion, a relaxation run assuming no sliding and a constant ice temperature of 10 C is performed to avoid spurious noise that arises from errors and biases in the data sets. To ensure that the relaxed geometry does not deviate too much from the observed geometry, the relaxation is conducted over 1 year. However, while inverse modelling is well established for estimating basal properties, the temperature field is difficult to constrain without performing an interglacial thermal spin-up. Therefore, we rely on a temperature field that was obtained by a hybrid approach between palaeoclimatic thermal spin-up and basal-friction inversion. This method was developed for the AWI contribution in initMIP-Greenland (Goelzer et al.2018) and further improved in Rückamp et al. (2018) by using the geothermal flux pattern from Greve (2005, scenario hf pmod2). Here, we initialize the ice rheology on the four employed G4000–G750 grids by interpolating the 3-D temperature and water content fields from the hybrid spin-up in Rückamp et al. (2018). The basal-melting rates of grounded ice are equivalently interpolated. During all transient runs, we neglect an evolution of the thermal field. This is justified as it was shown by Seroussi et al. (2013) and Goelzer et al. (2018, see submissions AWI-ISSM1 and 2) that the temperature field and its change have a negligible effect on century timescale projections of the GrIS.

The main ingredient of the initialization is the inversion to infer the basal-friction coefficient k2 in Eq. (1). This approach minimizes a cost function that measures the misfit between observed and modelled horizontal velocities (Morlighem et al.2010). The cost function is composed of two terms which fit the velocities in fast- and slow-moving areas. A third term is a Tikhonov regularization to avoid oscillations. The parameters for weighting the three contributions to the cost function are taken from Seroussi et al. (2013). We leverage horizontal surface velocities from the MEaSUREs project (Joughin et al.2016, 2018), as the data set with almost no gaps over the GrIS is suitable for basal-friction inversion.

The assigned reference year is 1990. This date is not in agreement with the timestamps of the BedMachine data set (reference time is 2007) and the MEaSUREs velocity data set (temporal coverage from 2014 to 2018). However, we ignore the contemporaneity requirement in the inversion approach and put more weight on starting the projections at the end of the assumed GrIS steady-state period (e.g. Ettema et al.2009). All transient simulations start from the initial state; that means, we do not perform a subsequent relaxation run to bring the model to a steady state (see Sect. 2.6).

2.4 Historical and control experiments

In both control experiments (ctrl and ctrl_proj), the SMB and ice sheet mask remain unchanged from the reference year according to the ISMIP6 protocol. To advance the model from the reference time to the projection start date, the historical scenario is needed. During the historical period, yearly cumulative SMB is taken from the RACMO2.3p2 product (Noël et al.2018) for the years from 1990 to 2015. For simplicity, the ice sheet extent remains unchanged from the reference year. This is a crude approach, but representing the historical mass loss accurately was not a strong priority for our experimental setup. As the ice front is not moving in these three scenarios, ice discharge Q equals calving D.

2.5 Future forcing experiments

It is beyond the scope of this paper to present the details of the ISMIP6 protocol and experimental design. Therefore, we aim to briefly outline the external-forcing approach. Further details are given in Goelzer et al. (2020a), Nowicki et al. (2020a), Fettweis et al. (2020) and Slater et al. (2019, 2020).

Table 2Summary of projection experiments based on MIROC5 RCP 8.5 climate data.

Download Print Version | Download XLSX

As we aim to study the effect of grid resolution on ice mass changes, we run the future scenarios based on climate data from one single GCM. The GCM MIROC5 was selected as it performs well in the historical period and represents a plausible climate near Greenland (Barthel et al.2020). The GCM output is used to separately derive ISM forcing for the interaction with the atmosphere and the ocean. We set up experiments where both external forcings are considered; these scenarios are termed “full” in the following (RCP8.5-Rlow, RCP8.5-Rmed, RCP8.5-Rhigh). In addition, we perform simulations where either the atmospheric forcing (RCP8.5-Rnone) or the marine-terminating outlet glacier retreat (OO-Rmed, OO-Rhigh) is at play. The conducted projection experiments and the corresponding experiment labels used in this study are summarized in Table 2 and are explained in the following sections.

2.5.1 Atmospheric forcing

ISMIP6 provides surface-forcing data sets for the GrIS based on CMIP GCM simulations. The GCM output is dynamically downscaled through the higher-resolution regional climate model (RCM) MAR v3.9 (Fettweis et al.2017). The latter allows for the capturing of narrow regions at the periphery of the Greenland ice sheet with large SMB gradients, which are likely not captured by the GCMs. The climatic SMB that is used as future-climate forcing reads

(2) SMB clim ( x , y , t ) = SMB ref ( x , y ) + Δ SMB ( x , y , t ) + SMB dyn ( x , y , t ) ,

with the anomaly defined as

(3) Δ SMB ( x , y , t ) = SMB ( x , y , t ) GCM–MAR - SMB ( x , y ) GCM–MAR 1960 - 1989 ,

where SMB(x,y,t)GCM–MAR is the direct output of MAR using the GCM climate data and SMB(x,y)GCM–MAR1960-1989 the corresponding mean value over the reference period (from January 1960 to December 1989). As the reference SMB field SMBref(x,y), we choose the downscaled RACMO2.3p2 product (Noël et al.2018) whereby a model output was averaged for the period 1960–1990. This period is chosen as the ice sheet is assumed to be close to a steady state in this period. (e.g. Ettema et al.2009). The SMB deduced by MAR is processed on a fixed topography (off-line); consequently local climate feedback processes due to the evolving surface in the projection experiments are not captured. The SMB height–elevation feedback is considered with a dynamic correction SMBdyn to the SMBclim following Franco et al. (2012):

(4) SMB dyn ( x , y , t ) = dSMBdz ( x , y , t ) × ( z s ( x , y , t ) - z ref ( x , y ) ) .

The surface elevation changes are taken from the ISM elevation zs(x,y,t) while running the simulation and the corresponding ISM reference elevation zref(x,y) from the initialization state. The yearly patterns of ΔSMB(x,y,t) and dSMBdz(x,y,t) are provided by ISMIP6. A cumulative SMB anomaly over the projection period is shown in Fig. 4a.

Figure 4Atmospheric and oceanic forcing. (a) Spatial pattern of the cumulative (2015–2100) SMB anomaly based on MIROC5 RCP 8.5 and downscaled with MAR (Fettweis et al.2020). (b) Retreat of marine-terminating outlet glaciers in the northwestern region under RCP8.5-Rhigh scenario. Purple areas indicate retreated areas in 2100. Region subsets are shown in Fig. 2.

2.5.2 Oceanic forcing

For the oceanic forcing we rely on the empirically derived outlet glacier parameterization retreat by Slater et al. (2019, 2020). This method circumvents the problem of employing a physically based calving law and frontal-melting rates based on GCM output. When employing this parameterization to the calving front, retreat and advance of marine-terminating outlet glaciers is directly prescribed as a yearly series of ice front positions. (i.e. it is not a result of the ice velocity at the calving front, calving rate and frontal melt that are used to simulate the calving-front position). Here, the binary retreat masks (i.e. ice-covered and non-ice-covered cells) are interpolated to the native grid by nearest-neighbour interpolation. Retreat occurs once a cell is fully emptied.

Though this parameterization is a strong simplification, it builds on projected sub-marine melting taking into account changes in ocean temperature and surface meltwater runoff from a GCM. The parameterization is not applied to the individual glaciers but over a predefined geographical region. Based on the numerous retreat trajectories, a medium retreat scenario as the trajectory with the median retreat in 2100 is defined. To cover uncertainty by this approach, low- and high-retreat scenarios are defined as the trajectories with the 25th and 75th percentile retreats in 2100. In the following, these retreat scenarios are termed Rlow, Rmed and Rhigh (Table 2). The retreat mask for RCP8.5-Rhigh in 2100 is exemplarily shown in Fig. 4b. Please note that the future-projection experiment RCP8.5-Rnone experiences no retreat of marine-terminating outlet glaciers.

2.6 Comparability of experiments

Central questions about resolution dependence are always as follows: how comparable are the results, and what is controlling the results? The presented initialization procedure and involved parameters are achieved for the high-resolution simulations (G750). The simulations with a coarser resolution would probably require other parameters, e.g. to obtain a better result for observational targets or to achieve a reduced model drift. However, we decided here to keep model parameters (e.g. inversion parameters) and parameterizations (e.g. subgrid scheme at grounding line) unchanged for all simulations. Similarly to the retreat masks, we rely on binary retreat masks for all adopted resolutions although the ISMIP6 protocol requests a subgrid scheme for coarse-resolution models. On the one hand this strategy simply assumes that the results are comparable as they build on the same basis. On the other hand it avoids exploring parameter spaces which are out of the scope of this study.

For the geometric input we are following the same strategy. It is always a compromise between matching the observed geometry or being closer to a steady state. Here, we put more weight on having the initial geometry closer to the observed geometry. Therefore, we directly started the historical run after the inversion, and no further relaxation run is performed to bring the model to a steady state. As the model is likely not in a steady state at the initial state, we expect a model drift in the transient runs which would not be the case for models that perform a relaxation towards a steady state after the inversion.

3 Results

3.1 Historical scenario

To evaluate the modelling decisions pertaining to the initialization, the state of the ice sheet at the end of the historical period is compared to observations. Due to the sparseness and limited temporal and spatial coverage of available observations, we rely on the BedMachine v3 (150 m grid spacing) and MEaSUREs (250 m grid spacing) data sets for ice thickness and surface velocity, respectively. As these data are used in the data assimilation and inversion, velocity and thickness are not independent quantities. However, during the historical period the ice sheet state is altered by the boundary conditions and external forcing. Therefore, the following evaluation attempts to quantify differences from the model configurations on the ISMIP6 projection start date.

Since the results are qualitatively similar for each grid simulation (Figs. S1, S2 and S3 in the Supplement), the surface velocity field of the G750 simulation is exemplarily shown in Fig. 5a. A consequence of the employed basal-friction inversion is the high fidelity in simulating the observed-velocity field indicated by a low root mean square error (RMSE; Fig. 5b). Notable is the decreasing RMSE with increasing spatial resolution. At the end of the historical experiment the RMSE is increased compared to the initialization due to geometric and velocity adjustments over the course of the experiment. However, the ice-sheet-wide RMSE of each model version is very similar, but in the areas of fast-flowing outlet glaciers (observed velocity > 200 m a−1) differences are more evident: the G4000 and G750 simulations yield RMSE=150ma-1 and RMSE=80ma-1, respectively. Note that these values are not identical to those given in Goelzer et al. (2020a), as the evaluation here is based on a different subsampling method. A mean signed difference (MSD) reflects a stronger underestimation of the simulated velocities with coarser resolution. The underestimation of prominent outlet glaciers for the G4000 setup is demonstrated in the spatial pattern of velocity differences (Fig. S4 in the Supplement). With increasing resolution, the difference pattern becomes more heterogeneous. Although barely visible, the G750 setup provides an interesting signature at narrowly confined outlet glaciers: generally, the velocities in the main trunk are underestimated, while beneath the shear margin velocities are overestimated. This might be due to the fact, that the employed resolution is not able to resolve the sharp velocity jump across the shear margin.

Figure 5Simulation results and error estimate of model output at the end of the historical run compared to observations. (a) Simulated surface velocity of the GrIS (m a−1) from the G750 simulation. The grey silhouette shows the Greenland land mask from BedMachine v3. Thin black lines show the grounding line. (b) Root mean square error (RMSE) of the horizontal velocity magnitude compared to MEaSUREs. (c) RMSE of ice thickness (not corrected with ctrl) compared to BedMachine v3. In panels (b) and (c), light and dark colours represent diagnostics at the initialization and end of the historical experiment, respectively. The diagnostics have been calculated on the regular MEaSUREs and BedMachine grids, respectively.

A similar evaluation for the thickness is performed. The ice-sheet-wide RMSE of ice thickness depicts the qualitatively similar grid-dependent behaviour as the velocity evaluation (Fig. 5c). Similarly, the RMSE shows larger differences in the fast-flow regions: the G4000 and G750 simulations yield RMSE=126 m and RMSE=45 m, respectively. In this region, the MSD indicates underestimation of ice thicknesses with coarser resolution. Spatial patterns of the thickness differences over the course of the historical experiment are shown in Fig. S5 in the Supplement.

The stored volumes, ice extent and spatially integrated SMB are among all grid simulations rather similar (V=7.28mSLE±0.2%, A= 1.787 × 106km2± 0.7 %, SMB = 375 Gt a−1± 0.2 %). However, the underestimation of velocities and ice thicknesses in the coarser-resolution models is confirmed by the temporal mean of the ice discharge in the historical period. The intrinsically simulated ice discharge Q yields 207 to 341 Gt a−1 for the G4000 and G750 simulations, respectively.

3.2 Sea-level contribution

In the following the transient effect of spatial resolution on ice volume evolution for the future-climate experiments is studied. The change in ice mass loss is expressed as sea-level contributions. Therefore, the simulated volume above flotation is converted into the total amount of global sea-level equivalent by assuming a constant ocean area of 3.618 × 108km2. In the following, the mass losses in the projection experiments are corrected with the ctrl run with respect to the reference time. For all conducted projection experiments, the determined GrIS mass losses as a function of time are shown in Fig. 6 and listed in Table 3.

Figure 6Projected sea-level contribution of the Greenland ice sheet based on MIROC5 RCP 8.5 climate data (a). Coloured lines indicate the different grid resolutions employed, while the individual scenarios are indicated with different line styles. The mass loss trends are corrected with the ctrl run relative to the reference time. The grey-shaded box shows the historical period. (b) Enlargement of the RCP8.5-Rlow, RCP8.5-Rmed, RCP8.5-Rhigh and RCP8.5-Rnone scenarios. (c) Enlargement of the OO-Rmed and OO-Rhigh scenarios.


Table 3Modelled mass change (mm SLE) in future experiments for all experiments.

 Numbers for G1000 and G750 are different compared to Goelzer et al. (2020a) as they are differently calculated (e.g. considered ocean area, native vs. interpolated
grid resolution).

Download Print Version | Download XLSX

As we have not initialized our model to be at a steady state, the transient response in the ctrl experiment (thin coloured lines in Fig. 6) should not be interpreted as a prediction of actual future behaviour; the ctrl run rather confirms that each model has achieved a high degree of equilibration, which is reflected by a low rate of volume change. As the initialization states are presumably different across the employed grids, we expect a different response of the ice sheet as it is likely not in equilibrium with the applied SMB and ice flux divergence. The simulated ice mass evolution shows for all models a mass gain for the 111-year ctrl experiment ranging between 28 and 2 mm. With increasing resolution, the drift gets smaller and is minimal for the G1000 and G750 simulations. Although projections are corrected with the ctrl run, the higher drift needs caution when interpreting the results as it has, e.g. a consequence on the SMB height–elevation feedback. The higher mass gain rates of the coarser resolutions in the ctrl simulation are due to the lower ice discharge rates (Sect. 3.1). Although the integrated signal in ice mass change is generally small, the spatial patterns reveal an ice thickness imbalance of up to hundreds of metres over the ctrl period (Fig. 7). Imposing an SMB correction to suppress the thickness imbalance would be feasible for maintaining a small drift. However, this is avoided here to enable a clean comparison between the four model versions and to leave the ice dynamics some degree of freedom. Moreover, the mass trends represent an important diagnostic. Comparing the ice thickness changes reveal distinct differences between the grid-resolution simulations (Fig. 7). For example, at the end of the ctrl run at some western and northwestern locations at the margin, the G4000 simulation exhibits thickening while the G750 reveals thinning. Another example is simulated at the southwestern margin, where extensive thickening prevails in all simulations but reaches farther inland in the coarser resolutions. However, from these figures, it becomes clear that positive and negative thickness changes partially compensate for each other, resulting in a low model drift.

Figure 7Difference in ice thickness between 2100 and 2015 for the ctrl run. (a) G4000 simulation and (b) G750 simulation. The grey silhouette shows the Greenland land mask from BedMachine v3. Positive values represent thickening, and negative values show thinning. Thin yellow line shows the grounding line in the year 2100.

Depending on the projection scenario, the GrIS will lose ice corresponding to an SLE of between 19 mm (or 108 mm excluding OO-Rmed and OO-Rhigh) and 137 mm. For the future-climate scenarios including atmospheric forcing a gradual increase in mass loss until the end of this century is simulated, indicating accelerating mass loss for a high-emission scenario. For RCP8.5-Rmed the mass loss reaches about 125.3 mm in 2100 (mean over G4000, G2000, G1000 and G750 results). The uncertainty quantification in the oceanic forcing results in a mean sea-level contribution that is 7.1 % less and 5.4 % greater for the RCP8.5-Rlow and RCP8.5-Rhigh scenarios, respectively. When no calving-front retreat is at play, i.e. the RCP8.5-Rnone scenario, the projected mean mass loss is approx. 105.0 mm, i.e  20 mm less compared to RCP8.5-Rmed. In contrast, the mean mass loss is considerably reduced to 26 and 37 mm in the OO-Rmed and OO-Rhigh experiment, respectively. Interestingly, a linear superposition of RCP8.5-Rnone and OO-Rmed leads to an overestimated mass loss of about 4.1 % for G4000 and 5.3 % for G750 compared to RCP8.5-Rmed where both external forcings are simultaneously at play; a linear superposition of RCP8.5-Rnone and OO-Rhigh leads to 4.5 % and 5.8 % overestimation. This is in line with earlier studies where this effect was already reported (Goelzer et al.2013; Fürst et al.2015)

Among all future projections a resolution-dependent impact on sea-level contribution is generally small compared to the total signal for our grids. In 2100, the spread in sea-level contribution is 6.4 mm in RCP8.5-Rhigh, 4.1 mm in RCP8.5-Rmed, 1.5 mm in RCP8.5-Rlow and 5 mm in RCP8.5-Rnone. Merely the OO-Rlow and OO-Rmed scenarios exhibit a spread of 10.7 mm and 13.6 mm, respectively, which is on the order of the absolute magnitude. A notable feature for all conducted simulations is that the sea-level contribution in each individual experiment converges with increasing resolution.

Figure 8 summarizes the qualitative behaviour of each experiment as a function of grid resolution. Note that the sea-level contribution in each experiment is normalized by its maximum. The finer resolutions tend to produce more mass loss in 2100 for the RCP8.5-Rmed, RCP8.5-Rhigh, OO-Rmed and OO-Rhigh experiments. An inverse behaviour is determined for the RCP8.5-Rnone experiment. The trend in the RCP8.5-Rlow experiment is not clear. The RCP8.5-Rnone, OO-Rmed and OO-Rhigh experiments unveil a linear behaviour as a function of grid size with regression slopes of m=1.50mmkm-1, m=-3.27mmkm-1 and m=-4.18mmkm-1, respectively. The trend in the full RCP8.5-Rlow, RCP8.5-Rmed and RCP8.5-Rhigh scenarios is not consistent: RCP8.5-Rmed and RCP8.5-Rhigh show a peak in mass loss at the finest resolution, whereas a peak in mass loss is attained in the G2000 simulation for RCP8.5-Rlow. For the latter, it is worth mentioning that the variations across the different grid simulations are less than 1.2 %. However, an intriguing effect of the conducted simulations remains the opposite behaviour of the RCP8.5-Rnone and e.g. RCP8.5-Rhigh scenarios. In the following section, we study this effect by analysing the mass partition to obtain a more in-depth insight into the role of atmospheric and oceanic forcing on grid resolution. It is worth mentioning that the qualitative behaviour of the detected grid-dependent mass loss remains similar when the projections are corrected with the ctrl_proj experiment (Table 3).

Figure 8Projected sea-level contribution in 2100 of the Greenland ice sheet as a function of the horizontal grid size. Values are normalized to the maximum of each experiment (coloured lines). Note the logarithmic scale of the x axis.


3.3 Mass partitioning

The relative mass loss partitioning in 2100 is shown in Fig. 9 to explore the role of the grid resolution in each experiment. The bars indicate the relative importance to sea-level contribution of ice dynamic changes in the projections. The dynamic contribution (composition of front retreat and ice discharge) is calculated as the residual of the total mass change and the integrated SMB anomaly. The remainder explains the part of SMB. The overall picture reveals that for experiments that include the atmospheric forcing the SMB anomaly is the governing forcing regardless of the grid resolutions. However, the importance of the dynamic contribution increases with larger prescribed retreat rates of outlet glaciers; i.e. G750 with RCP8.5-Rhigh on the upper end shows the highest importance of dynamic contribution with up to  28.4 %. On the lower end, the RCP8.5-Rnone shows diminished importance of dynamic contribution (< 5 %). In the OO-Rmed and OO-Rhigh scenarios, the mass loss is dominated by dynamic contribution. Concerning the grid resolution, the importance is on an equal level and exceeds 100 %. The negative importance of SMB stems from the fact that the glacier retreat is cutting off regions at the ice sheet margin where the static SMB is low.

Figure 9Mass loss partitioning for the conducted experiments. The bars indicate the relative dynamic contribution to sea level, calculated as the residual of the total of the mass change and the integrated SMB anomaly. The residual is a composition of front retreat and ice discharge.


In the full experiments RCP8.5-Rlow, RCP8.5-Rmed and RCP8.5-Rhigh, an increase in resolution enhances the importance of the dynamic contribution. For the G750 simulation it is  3 %, 5 % and 6 % higher for RCP8.5-Rlow, RCP8.5-Rmed and RCP8.5-Rhigh, respectively, compared to G4000. Curiously, the opposite behaviour is observed for the RCP8.5-Rnone experiment, where a finer resolution damps the importance of the dynamic contribution; G4000 yields a 4.9 % dynamic contribution, whereas G750 yields a 2.9 % dynamic contribution.

The simulated inverse grid-resolution responses raise the question of the driving causes. Overall the time series of the SMB show a decline and only minor differences among the grid resolutions (Fig. 10a). At the end of the projection, the cumulative SMB is 2.1 % and 2.6 % lower in the G4000 simulation for RCP8.5-Rnone and RCP8.5-Rhigh, respectively, compared to G750. These differences could be explained by different evolution of ablation areas at the margin and the SMB height–elevation feedback, in particular affected by the ctrl run, among all grid-resolution setups. In contrast, the cumulative ice discharge for these settings reveals an opposing response in the RCP8.5-Rnone and RCP8.5-Rhigh scenarios and more relative differences between the grid resolutions (Fig. 10b and c). At least for G2000, G1000 and G750, the ice discharge in the RCP8.5-Rnone experiment decreases over the century; the decrease in G4000 is offset by a few decades and exhibits an increase early in the century. These reductions explain the grid dependence of the dynamic contribution as listed in the previous paragraph (RCP8.5-Rnone in Fig. 9). For RCP8.5-Rhigh, the ice discharge shows an increase consistently but is more enhanced in the finer resolutions. This finding corroborates the grid-dependent increase in the relative ice discharge importance (RCP8.5-Rhigh in Fig. 9). As the opposing differences in RCP8.5-Rnone and RCP8.5-Rhigh are prevailing in ice discharge, it can be concluded that resolving ice discharge on the different grids is a decisive factor here. The involved feedback are further explored by focusing on particular outlet glaciers in the next section.

Figure 10Time series of cumulative SMB anomaly and cumulative ice discharge Q. Colour scheme is as in previous figures. (a) Cumulative SMB anomaly for RCP8.5-Rnone. (b) Cumulative ice discharge Q for RCP8.5-Rnone. (c) Ice discharge for RCP8.5-Rhigh (solid lines) and RCP8.5-Rlow (dashed lines). Cumulative SMB anomaly for RCP8.5-Rhigh and RCP8.5-Rlow is qualitatively similar to RCP8.5-Rnone. Ice discharge is not corrected with the ctrl run.


3.4 Outlet glacier response

The fact that the centennial mass loss for the full experiments increases as the grid size reduces raises the question whether this is caused by ice dynamics alone, dominant feedback with surface mass balance or retreat, or other non-obvious factors. We conduct an in-depth analysis of numerous prominent outlet glaciers at the GrIS (Fig. S6 and table with analysis provided as separate SI). The responses of most of the outlet glaciers reveal the deduced grid-dependent behaviour where higher resolutions cause an enhanced discharge. This is exemplarily illustrated in Fig. 11a for Helheim Glacier. However, this behaviour is not visible at all selected outlet glaciers. The presented example demonstrates that the bedrock topography deviates significantly among the different grid resolutions. Generally, the bedrock topography of the coarser resolution is located above the bed of the finer resolution. This topographic effect is restricted to narrow confined outlet glaciers that obey a characteristic width on the order of a few kilometres. Outlet glaciers that have a larger characteristic width, such as Humboldt Glacier, reveal in our setups a comparable bedrock topography. These glaciers seem to have a qualitatively equal behaviour for glacier speed-up and change in ice discharge for all employed grid resolutions (Fig. 11b). This analysis demonstrates that adjacent glaciers that experience similar environmental conditions may behave differently because ice discharge is strongly controlled by glacier geometry.

Figure 11Response of outlet glaciers. Colour scheme is the same as in previous figures, and light to dark colours indicate the years 2015, 2070 and 2100. (a) Helheim Glacier under RCP8.5-Rhigh forcing. (b) Humboldt Glacier under RCP8.5-Rhigh forcing. (c) Store Glacier under RCP8.5-Rhigh forcing. (d) Store Glacier under RCP8.5-Rnone forcing. Upper rows show the transient behaviour of the ice discharge q, the middle rows the surface velocity v and lower rows the evolution of the ice geometry. In the lower rows, the grey-shaded area depicts the bedrock topography from the G750 simulation. The grey lines from dark to light indicate the bedrock topography from the G1000, G2000 and G4000 simulation. None of the quantities are corrected with the ctrl run. Distance is relative to an arbitrary point.


Glaciers that are converted from a marine-terminating to a land-terminating glacier by retreating out of the water form their own class. These glaciers are no longer subject to the retreat and show a collapse in ice discharge regardless of the grid resolution as illustrated for Store Glacier in Fig. 11c. The qualitative behaviour of the retreat seems to be similar as reported in Aschwanden et al. (2019, Fig. 4b therein), but the timing of the retreat is different. In our study, Store Glacier is unstable and retreats within this century out of the water, while in Aschwanden et al. (2019) Store Glacier is in a very stable position; the quick retreat sets in far beyond 2100 once the glacier loses contact with the bedrock high. This different response is related to the retreat parameterization employed that lacks information on the bedrock topography, such as topographic highs and lows.

RCP8.5-Rnone shows a distinct slowdown in ice velocities as illustrated in Fig. 11d for Store Glacier. Visible is a larger slowdown of the higher resolutions; the same behaviour holds for the ice discharge q. This is in line with the finding above, i.e. that the scenario RCP8.5-Rnone reveals reduced ice discharge (Fig. 10b).

4 Discussion

The simulations presented here show that the projected sea-level contribution is sensitive to the spatial resolution. The sensitivity effect depends on the climate forcing, with oceanic and atmospheric forcings showing opposite and non-trivial responses. The simulations have shown that the ice discharge to the ocean is a decisive factor here controlling the grid-dependent spread. As shown above, outlet glaciers respond differently to external forcing, depending on the employed grids and geometrical setting. In such a non-linear system examining a driving mechanism remains challenging. However, despite the somewhat heterogeneous response of outlet glaciers, the different scenarios tend to produce an overall trend in characteristic fields that explains the different responses.

The different responses in the full scenarios could be attributed to the ability to resolve bedrock topography and the interaction with basal sliding. Figure 12 illustrates spatial changes in the effective pressure and basal-sliding velocity for RCP8.5-Rhigh. A common characteristic for G750 is a stronger decrease in effective pressure, which is concentrated in areas where the finer grid shows a deeper bed of the marine portions compared to G4000. Due to the linear dependence of τb on effective pressure (Eq. 1), basal-sliding velocities increase more strongly in the finer resolution. This feedback is enhanced as the SMB perturbations lead to a decrease in ice thickness, hence to a decrease in the effective pressure. The transient evolution reveals further that thinning and acceleration propagate faster and farther upstream in the finer resolution. The higher signal propagation rates may have additional consequences on longer timescales as the surface melt is amplified by the positive surface mass balance–elevation feedback exposing the ice surface to higher air temperatures.

Figure 12(a, b) Difference in effective pressure between 2100 and 2015 for the southeast region. (c, d) Difference in basal velocity between 2100 and 2015 for the southeast region. Region subsets are shown in Fig. 2. Dark grey line indicates the initial ice extent. Thin black line indicating the grounding line is not visible as it falls together with the calving front. The grey silhouette shows the Greenland land mask from BedMachine v3.

It remains questionable if the widespread glacier acceleration is induced by the frontal stress perturbation instead of the decrease in the effective pressure. To isolate this effect we conduct a RCP8.5-Rhigh simulation (not shown) where the effective pressure is held constant to the historical level. This setup reveals a very limited acceleration of a few glaciers in the G4000 simulation; some show no response or even a slowdown. In the corresponding G750 simulation most of the outlet glaciers show a speed-up, but this effect is very localized and does not reach far upstream. Therefore, we conclude, that the pronounced decrease in the effective pressure along with the acceleration of outlet glaciers is a dominant mechanism controlling the grid-dependent spread.

In order to investigate whether the response behaviour is an effect of purely reducing the grid size, we repeated the OO-Rhigh and RCP8.5-Rhigh experiments with a G1000 simulation using regridded bedrock topography and a friction coefficient from the G4000 initial state (simulations are not shown). This setup adopts a high-resolution grid but omits detailed information from the high-resolution input data. Projected mass loss by this setup is closer to the G4000 simulation. It, therefore, demonstrates that a high model resolution alone is insufficient to explain the grid-dependent sea-level contribution. As a consequence, a driving cause of the grid-dependent behaviour arises from additional information in the input data. Therefore, we conclude that the grid-dependent behaviour is highly connected to the bedrock topography because the different models represent the bedrock topography quite differently.

Compared to the full and ocean-only scenario, the atmospheric-only scenario unveils a stronger mass loss for the coarse resolution. To some extent, this could be explained by a slightly lower SMB in the coarser resolution. However, the finer resolution produces a stronger reduction in ice discharge over the course of the experiment. Although for many of the outlet glaciers the effective pressure decreases (not as strongly as for the scenarios with considered retreat), there is instead a slowdown of most glaciers (Fig. 13c and d). Again, these differences are concentrated in areas where the finer grid shows deeper troughs. Curiously, the finer resolution is better able to resolve these details, but the velocity evolution causes an extra reduction in ice discharge compared to in the coarser resolution. This non-trivial response is illustrated with the change in driving stress, approximated as τd=ϱigh|gradzs| (Fig. 13a and b). Compared to 2015, the driving stress has locally decreased more in the finer resolution in 2100 compared to the coarser resolution; away from the marginal region, the driving-stress changes are on a comparable magnitude for all employed grids. These results indicate that the reduction in the effective pressure is outperformed by geometric adjustments in the RCP8.5-Rnone scenario. This experiment intended to omit an interaction of the glacier with a changing ocean forcing, but the assumption of a fixed calving front hinders outlet glaciers from adjusting freely to topographic changes. They, therefore, do not experience reduced buttressing or frontal stress perturbations which are necessary mechanisms to trigger widespread glacier acceleration (e.g. Bondzio et al.2017). In future studies, it might be desirable to allow the calving front to adjust although the oceanic forcing is held constant. Nevertheless, the simulations indicate that without a frontal stress perturbation an ensuing speed-up of the outlet glacier is not initiated. This highlights the importance of capturing calving events, i.e tracking the ice front position in numerical models, most accurately.

Figure 13(a, b) Difference in driving stress between 2100 and 2015 for the northwestern region. (c, d) Difference in basal velocity between 2100 and 2015 for the northwestern region. Region subsets are shown in Fig. 2. Dark grey line indicates the initial ice extent. Thin black line indicating the grounding line is not visible as it falls together with the calving front. The grey silhouette shows the Greenland land mask from BedMachine v3.

The inverse grid-dependent behaviour of the RCP8.5-Rnone, OO-Rmed and OO-Rhigh scenarios has some implications when interpreting the mass loss of the ice sheet. The combined scenarios demonstrate that in a particular case, the sea-level contribution is maximized for an intermediate resolution. Depending on the horizontal grid resolution, the competing tendencies of SMB and ice discharge are differently resolved. This finding seems to corroborate results by Aschwanden et al. (2019), where an intermediate resolution reveals the largest sea-level contribution.

A convergence of the grid-dependent estimates of sea-level contribution emerges around REShigh 1 km. This value corroborates Aschwanden et al. (2016) for capturing outlet glacier behaviour indicating an upper limit for horizontal grid resolution. However, the converging behaviour should be treated with some caution. We cannot exclude the possibility that a model resolution finer than 750 m would lead to results that deviate from the convergence. On the one hand, the 150 m horizontal grid spacing of the BedMachine v3 data set is much finer than our finest resolution of 750 m. As the retreat parameterization is insensitive to bed undulations, resolving the outlet glacier cross sections is important for accurately modelling ice discharge. Since the glacier cross sections are reasonably well approximated in G750, we do not expect that resolving the geometry in a higher resolution would drastically alter ice discharge rates. On the other hand, there are indications that at a resolution of 750 m the HO solution is not fully converged (Pattyn et al.2008). Adopting a higher resolution could have implications for the ice flow and hence for the evolution of ice discharge. Likewise, in Aschwanden et al. (2019, Fig. S4 therein), the ice discharge is shown to increase as the mesh resolution is increased and seems to converge below a resolution of  1800 m. However, the finer resolutions of 450 and 600 m seem to produce again a somewhat lower ice discharge. This might indicate that the underlying processes are not fully converged and still causing changes in mass loss trends.

Our grid-dependent results under atmospheric-only forcing correlate with the finding in Goelzer et al. (2018, Fig. 1) and Exp. C2 in Greve and Herzfeld (2013, Fig. 7a and b therein). Interestingly, the causes for the same behaviour seem to have different origins. In Goelzer et al. (2018) the effect is likely due to an overestimated ablation area (see also Goelzer et al.2020b), whereas in our study the effect is attributed to the change in ice discharge of the ice sheet. The cause for the grid-dependent behaviour in Greve and Herzfeld (2013) is not specified further. Still, it is worth mentioning that they report a much better agreement of simulated to observed surface velocities by increasing the spatial resolution. Our experiments with the considered retreat of outlet glaciers could not be compared to the S1, M2 and R8 experiments in Greve and Herzfeld (2013). On the one hand, the external-forcing approach differs. On the other hand, a grid-dependent behaviour in Greve and Herzfeld (2013) is not clear (except for the enhanced sliding experiment S1, where the finer-resolution setups show a higher response). However, please note that the comparison to Aschwanden et al. (2019) and Greve and Herzfeld (2013) is limited as the studies employ different flow approximations and treat basal friction differently. The underlying physics in each ISM probably depends differently on the resolution.

Besides the fixed calving front in the atmospheric-only scenario, further limitations of our study must be noted. The spread of projected sea-level contribution among all grids is likely subject to the chosen type of friction law. The choice of the basal-friction law used in ISMs remains a matter of debate (Stearns and van der Veen2018; Minchew et al.2019) and a potential source of uncertainty in sea-level projections (e.g. Brondex et al.2019). Compared to our used type of friction law, there are some indications that friction laws satisfying an upper bound for the basal drag are more reliable (e.g. Schoof2005; Gagliardini et al.2007; Leguy et al.2014; Joughin et al.2019). In brief, these types of friction laws invoke a switch for the friction regimes (low and high N, respectively) so that the influence of the effective pressure on the basal drag at slow ice flow is vanishing. It would be most interesting to evaluate their sensitivity to the horizontal resolution for projections on centennial timescales.

Another limitation concerns the choice of inversion parameters. We performed the inversions for basal parameters for each grid resolution individually but relying on the inversion parameters tuned for the high-resolution setups. Effectively, this results in an overall comparable pattern for the flow velocities (Figs. S1 and S2) and basal-friction coefficient (Fig. S3) for all grids. However, on smaller scales, the inversion approach produces significantly different k2 in many glacier basins. Recalling the relationship between N and k2, these different patterns are plausible but could potentially be a result from non-optimal inversion parameters. However, these different spatial patterns might make an additional contribution to the grid dependence of the simulations. In future studies, it will be worth investigating this influence only, e.g. by tuning the inversion parameters for each grid separately to find the optimal parameters.

However, the simulations conducted here reveal a grid-dependent spread in the full scenarios ranging between 1.2 % and 5.3 %, which is of a comparable magnitude to the surface mass balance–elevation feedback (Eq. 4). The latter is recognized as an important mechanism and accounts for an additional sea-level contribution of about 6 %–8 % (Goelzer et al.2020a). A feedback that we have not considered is the enhanced surface melt influencing the basal conditions and calving by filling up crevasses. Given that we greatly simplify the representation of N, the effect of a reorganization of basal conditions (in either way) or the effect of increased availability of water due to increasing surface melt on basal sliding is suppressed. To overcome this limitation, an adequate subglacial hydrology model could be invoked, even if not considering seasonality. A subglacial hydrology model has shown the localized effect on N, which likely has consequences for the spread between the employed grid resolutions (e.g. Werder et al.2013; de Fleurian et al.2016; Rathmann et al.2017; Sommers et al.2018; Beyer et al.2018; Neckel et al.2020).

A feedback that is not fully covered in our simulations is shear margin weakening and its influence on the evolution of flow velocities. Although the shear margins are weakly developed in the simulations (more pronounced in the finer resolutions), it is expected that a thermo-mechanical coupling could further weaken the shear margins as a response to frontal stress perturbations (e.g. Bondzio et al.2017). Such a coupling would increase the widespread inland flow acceleration and enhance the rate of mass loss. However, the change in τb and τd is very pronounced in and around the main trunks and quite different for the adopted grids (see Figs. S20 and S21 in the Supplement). These patterns exemplify the need for resolving the shear margins in a particularly high resolution, which we have not fully accomplished in this study as shear margins are becoming, in numerous cases, subgrid phenomena. This may be a reason for under-representing glacier velocities inside the main trunk and overestimating velocities outside the main flow as apparent in G750. This effect shall be addressed in further studies, in which ideally higher resolutions are employed or error estimators are engaged (e.g. dos Santos et al.2019).

5 Conclusions

We applied the three-dimensional finite-element higher-order model ISSM to the Greenland ice sheet to simulate the future response under climatic changes specified by the ISMIP6 protocol. The sensitivity of mass changes to the spatial resolution is tested by employing four different grids with a varying horizontal resolution ranging from 4 to 0.75 km at fast-flowing outlet glaciers. The simulations reveal up to  5.3 % more sea-level rise compared to the coarser resolution in the full scenario RCP8.5-Rhigh and  3.2 % for RCP8.5-Rmed. In scenarios where a change in SMB is omitted and only outlet glacier retreat is at play, the finer resolutions produce significantly more mass loss (up to 33 %). When no retreat is enforced, the sensitivity of the grid dependence exhibits an inverse behaviour; i.e. the coarser resolutions produce more mass loss. This finding is important to recognize for ice sheet models that have SMB as the dominant mass loss driver.

The results presented underline the importance of resolving the bedrock topography accurately. Areas with simple and low bedrock undulations experience a similar response in all model resolutions. In areas with complex and high bedrock undulations striking differences between the employed grids emerge. A mechanism that exerts an important control on the resolution-dependent spread is basal sliding predominantly in marine portions of outlet glaciers. Since we rely on a greatly simplified effective-pressure parameterization, further work is needed to prove the robustness of this conclusion.

Given the strong interaction of the bedrock topography with sliding, it is obvious that the major outlet glacier should be surveyed with the latest radar technology to obtain a substantially improved survey of the bedrock topography, the area of expected retreat and connected areas further upstream. This, in turn, requires ice sheet models ready to resolve these areas in terms of grids and physics adequately.

Code and data availability

Simulation results on the native grids described in this paper will be made publicly available with a digital object identifier: (Rückamp et al.2020). The forcing data sets are available through (ISMIP62015, last access: 21 August 2020). The ice flow model ISSM is open source and freely available at (last access: 21 August 2020; Larour et al.2012). Here ISSM version 4.16 is used.


The supplement related to this article is available online at:

Author contributions

MR conducted the study supported by the other authors. MR set up the ISSM model and ran the experiments. AH analysed the results for the individual glaciers. HG calculated the retreat masks. MR wrote the manuscript together with the other authors.

Competing interests

The authors declare that they have no conflict of interest.

Special issue statement

This article is part of the special issue “The Ice Sheet Model Intercomparison Project for CMIP6 (ISMIP6)”. It is not associated with a conference.


We thank the Climate and Cryosphere (CliC) effort, which provided support for ISMIP6 through sponsoring of workshops, hosting the ISMIP6 website and wiki, and promoting ISMIP6. We acknowledge the World Climate Research Programme, which, through its Working Group on Coupled Modelling, coordinated and promoted CMIP5 and CMIP6. We thank the climate-modelling groups for producing and making available their model output, the Earth System Grid Federation (ESGF) for archiving the CMIP data and providing access, the University at Buffalo for ISMIP6 data distribution and upload, and the multiple funding agencies who support CMIP5 and CMIP6 and ESGF. We thank the ISMIP6 steering committee, the ISMIP6 model selection group and ISMIP6 data set preparation group for their continuous engagement in defining ISMIP6 and for all their efforts.

We thank Thomas Kleiner (AWI) and Ralf Greve (ILTS) for continuous and fruitful discussions on simulations. We would like to thank Natalja Rakowsky, Malte Thoma and Bernadette Fritzsch for maintaining excellent computing facilities at AWI and DKRZ. The high-resolution simulations were run on the DKRZ HPC system Mistral under grant ab1073. We also acknowledge the outstanding support of the ISSM team. We thank the editor Robin Smith and the two anonymous reviewers for their helpful suggestions to improve the manuscript.

This is ISMIP6 contribution no. 17.

Financial support

Martin Rückamp has been supported by the Helmholtz Climate Initiative REKLIM (regional climate change). Heiko Goelzer has received funding from the programme of the Netherlands Earth System Science Centre (NESSC), financially supported by the Dutch Ministry of Education, Culture and Science (OCW; grant no. 024.002.001).

The article processing charges for this open-access
publication were covered by a Research
Centre of the Helmholtz Association.

Review statement

This paper was edited by Robin Smith and reviewed by two anonymous referees.


Aschwanden, A., Fahnestock, M. A., and Truffer, M.: Complex Greenland outlet glacier flow captured, Nat. Commun., 7, 10524,, 2016. a, b, c

Aschwanden, A., Fahnestock, M. A., Truffer, M., Brinkerhoff, D. J., Hock, R., Khroulev, C., Mottram, R., and Khan, S. A.: Contribution of the Greenland Ice Sheet to sea level over the next millennium, Science Advances, 5, 6,, 2019. a, b, c, d, e, f, g, h

Barthel, A., Agosta, C., Little, C. M., Hattermann, T., Jourdain, N. C., Goelzer, H., Nowicki, S., Seroussi, H., Straneo, F., and Bracegirdle, T. J.: CMIP5 model selection for ISMIP6 ice sheet model forcing: Greenland and Antarctica, The Cryosphere, 14, 855–879,, 2020. a

Beckmann, A. and Goosse, H.: A parameterization of ice shelf-ocean interaction for climate models, Ocean Model., 5, 157–170,, 2003. a

Beyer, S., Kleiner, T., Aizinger, V., Rückamp, M., and Humbert, A.: A confined–unconfined aquifer model for subglacial hydrology and its application to the Northeast Greenland Ice Stream, The Cryosphere, 12, 3931–3947,, 2018. a

Bindschadler, R. A., Nowicki, S., Abe-Ouchi, A., Aschwanden, A., Choi, H., Fastook, J., Granzow, G., Greve, R., Gutowski, G., Herzfeld, U. C., Jackson, C., Johnson, J., Khroulev, C., Levermann, A., Lipscomb, W. H., Martin, M. A., Morlighem, M., Parizek, B. R., Pollard, D., Price, S. F., Ren, D., Saito, F., Sato, T., Seddik, H., Seroussi, H., Takahashi, K., Walker, R., and Wang, W. L.: Ice-sheet model sensitivities to environmental forcing and their use in projecting future sea level (the SeaRISE project), J. Glaciol., 59, 195–224,, 2013. a

Blatter, H.: Velocity and stress fields in grounded glaciers: a simple algorithm for including deviatoric stress gradients, J. Glaciol., 41, 333–344,, 1995. a

Bondzio, J. H., Morlighem, M., Seroussi, H., Kleiner, T., Rückamp, M., Mouginot, J., Moon, T., Larour, E. Y., and Humbert, A.: The mechanisms behind Jakobshavn Isbræ's acceleration and mass loss: A 3-D thermomechanical model study, Geophys. Res. Lett., 44, 6252–6260,, 2017. a, b, c

Brondex, J., Gillet-Chaulet, F., and Gagliardini, O.: Sensitivity of centennial mass loss projections of the Amundsen basin to the friction law, The Cryosphere, 13, 177–195,, 2019. a

Budd, W. F., Jenssen, D., and Smith, I. N.: A Three-Dimensional Time-Dependent Model of the West Antarctic Ice Sheet, Ann. Glaciol., 5, 29–36,, 1984. a

Choi, Y., Morlighem, M., Rignot, E., Mouginot, J., and Wood, M.: Modeling the Response of Nioghalvfjerdsfjorden and Zachariae Isstrøm Glaciers, Greenland, to Ocean Forcing Over the Next Century, Geophys. Res. Lett., 44, 11071–11079,, 2017. a

Cogley, J., Hock, R., Rasmussen, L., Arendt, A., Bauder, A., Braithwaite, R., Jansson, P., Kaser, G., Möller, M., Nicholson, L., and Zemp, M.: Glossary of Glacier Mass Balance and Related Terms, Tech. Rep. 86, IHP-VII Technical Documents in Hydrology, UNESCO-IHP, Paris, France, 2011. a

de Fleurian, B., Morlighem, M., Seroussi, H., Rignot, E., van den Broeke, M. R., Kuipers Munneke, P., Mouginot, J., Smeets, P. C. J. P., and Tedstone, A. J.: A modeling study of the effect of runoff variability on the effective pressure beneath Russell Glacier, West Greenland, J. Geophys. Res.-Earth, 121, 1834–1848,, 2016. a

dos Santos, T. D., Morlighem, M., Seroussi, H., Devloo, P. R. B., and Simões, J. C.: Implementation and performance of adaptive mesh refinement in the Ice Sheet System Model (ISSM v4.14), Geosci. Model Dev., 12, 215–232,, 2019. a

Enderlin, E. M., Howat, I. M., Jeong, S., Noh, M.-J., Angelen, J. H., and Broeke, M. R.: An improved mass budget for the Greenland ice sheet, Geophys. Res. Lett., 41, 866–872,, 2014. a

Ettema, J., van den Broeke, M. R., van Meijgaard, E., van de Berg, W. J., Bamber, J. L., Box, J. E., and Bales, R. C.: Higher surface mass balance of the Greenland ice sheet revealed by high-resolution climate modeling, Geophys. Res. Lett., 36, 12,, 2009. a, b

Eyring, V., Bony, S., Meehl, G. A., Senior, C. A., Stevens, B., Stouffer, R. J., and Taylor, K. E.: Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization, Geosci. Model Dev., 9, 1937–1958,, 2016. a

Fettweis, X., Box, J. E., Agosta, C., Amory, C., Kittel, C., Lang, C., van As, D., Machguth, H., and Gallée, H.: Reconstructions of the 1900–2015 Greenland ice sheet surface mass balance using the regional climate MAR model, The Cryosphere, 11, 1015–1033,, 2017. a

Fettweis, X., Hofer, S., Krebs-Kanzow, U., Amory, C., Aoki, T., Berends, C. J., Born, A., Box, J. E., Delhasse, A., Fujita, K., Gierz, P., Goelzer, H., Hanna, E., Hashimoto, A., Huybrechts, P., Kapsch, M.-L., King, M. D., Kittel, C., Lang, C., Langen, P. L., Lenaerts, J. T. M., Liston, G. E., Lohmann, G., Mernild, S. H., Mikolajewicz, U., Modali, K., Mottram, R. H., Niwano, M., Noël, B., Ryan, J. C., Smith, A., Streffing, J., Tedesco, M., van de Berg, W. J., van den Broeke, M., van de Wal, R. S. W., van Kampenhout, L., Wilton, D., Wouters, B., Ziemen, F., and Zolles, T.: GrSMBMIP: Intercomparison of the modelled 1980–2012 surface mass balance over the Greenland Ice sheet, The Cryosphere Discuss.,, in review, 2020. a, b

Franco, B., Fettweis, X., Lang, C., and Erpicum, M.: Impact of spatial resolution on the modelling of the Greenland ice sheet surface mass balance between 1990–2010, using the regional climate model MAR, The Cryosphere, 6, 695–711,, 2012. a

Fürst, J. J., Goelzer, H., and Huybrechts, P.: Ice-dynamic projections of the Greenland ice sheet in response to atmospheric and oceanic warming, The Cryosphere, 9, 1039–1062,, 2015. a, b, c

Gagliardini, O., Cohen, D., Raback, P., and Zwinger, T.: Finite-element modeling of subglacial cavities and related friction law, J. Geophys. Res.-Earth, 112, F2,, 2007. a, b

Gillet-Chaulet, F., Gagliardini, O., Seddik, H., Nodet, M., Durand, G., Ritz, C., Zwinger, T., Greve, R., and Vaughan, D. G.: Greenland ice sheet contribution to sea-level rise from a new-generation ice-sheet model, The Cryosphere, 6, 1561–1576,, 2012. a, b

Gladstone, R., Schäfer, M., Zwinger, T., Gong, Y., Strozzi, T., Mottram, R., Boberg, F., and Moore, J. C.: Importance of basal processes in simulations of a surging Svalbard outlet glacier, The Cryosphere, 8, 1393–1405,, 2014. a

Glen, J. W.: The Creep of Polycrystalline Ice, P. Roy. Soc. Lond. A Mat., 228, 519–538,, 1955. a

Goelzer, H., Huybrechts, P., Fürst, J., Nick, F., Andersen, M., Edwards, T., Fettweis, X., Payne, A., and Shannon, S.: Sensitivity of Greenland Ice Sheet Projections to Model Formulations, J. Glaciol., 59, 733–749,, 2013. a, b, c

Goelzer, H., Nowicki, S., Edwards, T., Beckley, M., Abe-Ouchi, A., Aschwanden, A., Calov, R., Gagliardini, O., Gillet-Chaulet, F., Golledge, N. R., Gregory, J., Greve, R., Humbert, A., Huybrechts, P., Kennedy, J. H., Larour, E., Lipscomb, W. H., Le clec'h, S., Lee, V., Morlighem, M., Pattyn, F., Payne, A. J., Rodehacke, C., Rückamp, M., Saito, F., Schlegel, N., Seroussi, H., Shepherd, A., Sun, S., van de Wal, R., and Ziemen, F. A.: Design and results of the ice sheet model initialisation experiments initMIP-Greenland: an ISMIP6 intercomparison, The Cryosphere, 12, 1433–1460,, 2018. a, b, c, d, e, f, g, h

Goelzer, H., Nowicki, S., Payne, A., Larour, E., Seroussi, H., Lipscomb, W. H., Gregory, J., Abe-Ouchi, A., Shepherd, A., Simon, E., Agosta, C., Alexander, P., Aschwanden, A., Barthel, A., Calov, R., Chambers, C., Choi, Y., Cuzzone, J., Dumas, C., Edwards, T., Felikson, D., Fettweis, X., Golledge, N. R., Greve, R., Humbert, A., Huybrechts, P., Le clec'h, S., Lee, V., Leguy, G., Little, C., Lowry, D. P., Morlighem, M., Nias, I., Quiquet, A., Rückamp, M., Schlegel, N.-J., Slater, D. A., Smith, R. S., Straneo, F., Tarasov, L., van de Wal, R., and van den Broeke, M.: The future sea-level contribution of the Greenland ice sheet: a multi-model ensemble study of ISMIP6, The Cryosphere, 14, 3071–3096,, 2020a. a, b, c, d, e, f, g, h, i

Goelzer, H., Noël, B. P. Y., Edwards, T. L., Fettweis, X., Gregory, J. M., Lipscomb, W. H., van de Wal, R. S. W., and van den Broeke, M. R.: Remapping of Greenland ice sheet surface mass balance anomalies for large ensemble sea-level change projections, The Cryosphere, 14, 1747–1762,, 2020b. a

Greve, R.: Relation of measured basal temperatures and the spatial distribution of the geothermal heat flux for the Greenland ice sheet, Ann. Glaciol., 42, 424–432, 2005. a

Greve, R. and Herzfeld, U. C.: Resolution of ice streams and outlet glaciers in large-scale simulations of the Greenland ice sheet, Ann. Glaciol., 54, 209–220,, 2013. a, b, c, d, e, f, g, h

Iken, A.: The Effect of the Subglacial Water Pressure on the Sliding Velocity of a Glacier in an Idealized Numerical Model, J. Glaciol., 27, 407–421,, 1981. a

ISMIP6: ISMIP6 wiki page, available at: (last access 21 August 2020), 2015. a

Joughin, I., Smith, B., Howat, I., and Scambos, T.: MEaSUREs Multi-year Greenland Ice Sheet Velocity Mosaic, Version 1, NASA National Snow and Ice Data Center Distributed Active Archive Center, Boulder, Colorado, USA,, 2016. a, b

Joughin, I., Smith, B., and Howat, I.: A complete map of Greenland ice velocity derived from satellite data collected over 20 years, J. Glaciol., 64, 1–11,, 2018. a, b

Joughin, I., Smith, B. E., and Schoof, C. G.: Regularized Coulomb Friction Laws for Ice Sheet Sliding: Application to Pine Island Glacier, Antarctica, Geophys. Res. Lett., 46, 4764–4771,, 2019. a, b

Larour, E., Seroussi, H., Morlighem, M., and Rignot, E.: Continental scale, high order, high spatial resolution, ice sheet modeling using the Ice Sheet System Model (ISSM), J. Geophys. Res.-Earth, 117, F01022,, 2012 (ISSM source code available at: (last access: 21 August 2020). a, b, c, d

Leguy, G. R., Asay-Davis, X. S., and Lipscomb, W. H.: Parameterization of basal friction near grounding lines in a one-dimensional ice sheet model, The Cryosphere, 8, 1239–1259,, 2014. a, b

Lliboutry, L. and Duval, P.: Various isotropic and anisotropic ices found in glaciers and polar ice caps and their corresponding rheologies, Ann. Geophys., 3, 207–224, 1985. a

Minchew, B. M., Meyer, C. R., Pegler, S. S., Lipovsky, B. P., Rempel, A. W., Gudmundsson, G. H., and Iverson, N. R.: Comment on “Friction at the bed does not control fast glacier flow”, Science, 363, 6427,, 2019. a

Moon, T., Joughin, I., Smith, B., and Howat, I.: 21st-Century Evolution of Greenland Outlet Glacier Velocities, Science, 336, 576–578,, 2012. a

Morlighem, M., Rignot, E., Seroussi, H., Larour, E., Dhia, H. B., and Aubry, D.: Spatial patterns of basal drag inferred using control methods from a full-Stokes and simpler models for Pine Island Glacier, West Antarctica, Geophys. Res. Lett., 37, 14,, 2010. a, b, c

Morlighem, M., Williams, C. N., Rignot, E., An, L., Arndt, J. E., Bamber, J. L., Catania, G., Chauché, N., Dowdeswell, J. A., Dorschel, B., Fenty, I., Hogan, K., Howat, I., Hubbard, A., Jakobsson, M., Jordan, T. M., Kjeldsen, K. K., Millan, R., Mayer, L., Mouginot, J., Noël, B. P. Y., O'Cofaigh, C., Palmer, S., Rysgaard, S., Seroussi, H., Siegert, M. J., Slabon, P., Straneo, F., van den Broeke, M. R., Weinrebe, W., Wood, M., and Zinglersen, K. B.: BedMachine v3: Complete bed topography and ocean bathymetry mapping of Greenland from multibeam echo sounding combined with mass conservation, Geophys. Res. Lett., 44, 11051–11061,, 2017. a, b

Moss, R. H., Edmonds, J. A., Hibbard, K. A., Manning, M. R., Rose, S. K., van Vuuren, D. P., Carter, T. R., Emori, S., Kainuma, M., Kram, T., Meehl, G. A., Mitchell, J. F. B., Nakicenovic, N., Riahi, K., Smith, S. J., Stouffer, R. J., Thomson, A. M., Weyant, J. P., and Wilbanks, T. J.: The next generation of scenarios for climate change research and assessment, Nature, 463, 747–756,, 2010. a

Mouginot, J., Rignot, E., Bjørk, A. A., van den Broeke, M., Millan, R., Morlighem, M., Noël, B., Scheuchl, B., and Wood, M.: Forty-six years of Greenland Ice Sheet mass balance from 1972 to 2018, P. Natl. Acad. Sci. USA, 116, 9239–9244,, 2019. a

Neckel, N., Zeising, O., Steinhage, D., Helm, V., and Humbert, A.: Seasonal Observations at 79 N Glacier (Greenland) From Remote Sensing and in situ Measurements, Front. Earth Sci., 8, 142,, 2020. a

Nerem, R. S., Beckley, B. D., Fasullo, J. T., Hamlington, B. D., Masters, D., and Mitchum, G. T.: Climate-change-driven accelerated sea-level rise detected in the altimeter era, P. Natl. Acad. Sci. USA, 115, 2022–2025,, 2018. a

Noël, B., van de Berg, W. J., van Wessem, J. M., van Meijgaard, E., van As, D., Lenaerts, J. T. M., Lhermitte, S., Kuipers Munneke, P., Smeets, C. J. P. P., van Ulft, L. H., van de Wal, R. S. W., and van den Broeke, M. R.: Modelling the climate and surface mass balance of polar ice sheets using RACMO2 – Part 1: Greenland (1958–2016), The Cryosphere, 12, 811–831,, 2018. a, b

Nowicki, S., Bindschadler, R. A., Abe-Ouchi, A., Aschwanden, A., Bueler, E., Choi, H., Fastook, J., Granzow, G., Greve, R., Gutowski, G., Herzfeld, U., Jackson, C., Johnson, J., Khroulev, C., Larour, E., Levermann, A., Lipscomb, W. H., Martin, M. A., Morlighem, M., Parizek, B. R., Pollard, D., Price, S. F., Ren, D., Rignot, E., Saito, F., Sato, T., Seddik, H., Seroussi, H., Takahashi, K., Walker, R., and Wang, W. L.: Insights into spatial sensitivities of ice mass response to environmental change from the SeaRISE ice sheet modeling project II: Greenland, J. Geophys. Res.-Earth, 118, 1025–1044,, 2013. a

Nowicki, S., Goelzer, H., Seroussi, H., Payne, A. J., Lipscomb, W. H., Abe-Ouchi, A., Agosta, C., Alexander, P., Asay-Davis, X. S., Barthel, A., Bracegirdle, T. J., Cullather, R., Felikson, D., Fettweis, X., Gregory, J. M., Hattermann, T., Jourdain, N. C., Kuipers Munneke, P., Larour, E., Little, C. M., Morlighem, M., Nias, I., Shepherd, A., Simon, E., Slater, D., Smith, R. S., Straneo, F., Trusel, L. D., van den Broeke, M. R., and van de Wal, R.: Experimental protocol for sea level projections from ISMIP6 stand-alone ice sheet models, The Cryosphere, 14, 2331–2368,, 2020a. a, b, c

Nowicki, S., Payne, A. J., Abe-Ouchi, A., Agosta, C., Alexander, P., Albrecht, T., Asay-Davis, X., Barthel, A., Calov, R., Chambers, C., Choi, Y., Cullather, R., Cuzzone, J., Dumas, C., Edwards, T., Felikson, D., Fettweis, X., Goelzer, H., Gladstone, R., Golledge, N. R., Gregory, J. M., Greve, R., Hatterman, T., Hoffman, M. J., Humbert, A., Huybrechts, P., Jourdain, N. C., Kleiner, T., Larour, E., Le clec'h, S., Lee, V., Leguy, G., Lipscomb, W. H., Little, C. M., Lowry, D. P., Morlighem, M., Nias, I., Pattyn, F., Pelle, T., Price, S., Quiquet, A., Reese, R., Rückamp, M., Schlegel, N.-J., Seroussi, H., Shepherd, A., Simon, E., Slater, D., Smith, R., Straneo, F., Sun, S., Tarasov, L., Trusel, L. D., Breedam, J. V., van de Wal, R., van den Broeke, M., Winkelmann, R., Zhao, C., Zhang, T., and Zwinger, T.: Contrasting contributions to future sea level under CMIP5 and CMIP6 scenarios from the Greenland and Antarctic ice sheets, Geophys. Res. Lett., in review, 2020b. a

Nowicki, S. M. J., Payne, A., Larour, E., Seroussi, H., Goelzer, H., Lipscomb, W., Gregory, J., Abe-Ouchi, A., and Shepherd, A.: Ice Sheet Model Intercomparison Project (ISMIP6) contribution to CMIP6, Geosci. Model Dev., 9, 4521–4545,, 2016. a

Pattyn, F.: A new three-dimensional higher-order thermomechanical ice-sheet model: basic sensitivity, ice-stream development and ice flow across subglacial lakes, J. Geophys. Res.-Sol. Ea., 108, 2382,, 2003. a

Pattyn, F. and Durand, G.: Why marine ice sheet model predictions may diverge in estimating future sea level rise, Geophys. Res. Lett., 40, 4316–4320,, 2013. a

Pattyn, F., Perichon, L., Aschwanden, A., Breuer, B., de Smedt, B., Gagliardini, O., Gudmundsson, G. H., Hindmarsh, R. C. A., Hubbard, A., Johnson, J. V., Kleiner, T., Konovalov, Y., Martin, C., Payne, A. J., Pollard, D., Price, S., Rückamp, M., Saito, F., Souček, O., Sugiyama, S., and Zwinger, T.: Benchmark experiments for higher-order and full-Stokes ice sheet models (ISMIP–HOM), The Cryosphere, 2, 95–108,, 2008. a, b

Perego, M., Price, S., and Stadler, G.: Optimal initial conditions for coupling ice sheet models to Earth system models, J. Geophys. Res.-Earth, 119, 1894–1917,, 2014. a

Price, S. F., Payne, A. J., Howat, I. M., and Smith, B. E.: Committed sea-level rise for the next century from Greenland ice sheet dynamics during the past decade, P. Natl. Acad. Sci. USA, 108, 8978–8983,, 2011. a

Rathmann, N. M., Hvidberg, C. S., Solgaard, A. M., Grinsted, A., Gudmundsson, G. H., Langen, P. L., Nielsen, K. P., and Kusk, A.: Highly temporally resolved response to seasonal surface melt of the Zachariae and 79N outlet glaciers in northeast Greenland, Geophys. Res. Lett., 44, 9805–9814,, 2017. a

Rietbroek, R., Brunnabend, S.-E., Kusche, J., Schröter, J., and Dahle, C.: Revisiting the contemporary sea-level budget on global and regional scales, P. Natl. Acad. Sci. USA, 113, 1504–1509,, 2016. a

Rückamp, M., Falk, U., Frieler, K., Lange, S., and Humbert, A.: The effect of overshooting 1.5 C global warming on the mass loss of the Greenland ice sheet, Earth Syst. Dynam., 9, 1169–1189,, 2018. a, b, c

Rückamp, M., Greve, R., and Humbert, A.: Comparative simulations of the evolution of the Greenland ice sheet under simplified Paris Agreement scenarios with the models SICOPOLIS and ISSM, Polar Sci., 21, 14–25,, 2019a. a

Rückamp, M., Neckel, N., Berger, S., Helm, V., and Humbert, A.: Calving induced Speedup of Petermann Glacier, J. Geophys. Res.-Earth, 124, 216–228,, 2019b. a, b

Rückamp, M., Goelzer, H., and Humbert, A.: Results of “Sensitivity of Greenland ice sheet projections to spatial resolution in higher-order simulations: the AWI contribution to ISMIP6-Greenland using ISSM”, Zenodo,, 2020. a

Schoof, C.: The effect of cavitation on glacier sliding, P. Roy. Soc. A-Math. Phy., 461, 609–627,, 2005. a, b

Seroussi, H. and Morlighem, M.: Representation of basal melting at the grounding line in ice flow models, The Cryosphere, 12, 3085–3096,, 2018. a

Seroussi, H., Morlighem, M., Rignot, E., Khazendar, A., Larour, E., and Mouginot, J.: Dependence of century-scale projections of the Greenland ice sheet on its thermal regime, J. Glaciol., 59, 1024–1034,, 2013. a, b, c, d, e

Seroussi, H., Morlighem, M., Larour, E., Rignot, E., and Khazendar, A.: Hydrostatic grounding line parameterization in ice sheet models, The Cryosphere, 8, 2075–2087,, 2014. a

Shepherd, A., Ivins, E., Rignot, E., Smith, B., van den Broeke, M., Velicogna, I., Whitehouse, P., Briggs, K., Joughin, I., Krinner, G., Nowicki, S., Payne, T., Scambos, T., Schlegel, N., Geruo, A., Agosta, C., Ahlstrøm, A., Babonis, G., Barletta, V. R., Bjørk, A. A., Blazquez, A., Bonin, J., Colgan, W., Csatho, B., Cullather, R., Engdahl, M. E., Felikson, D., Fettweis, X., Forsberg, R., Hogg, A. E., Gallee, H., Gardner, A., Gilbert, L., Gourmelen, N., Groh, A., Gunter, B., Hanna, E., Harig, C., Helm, V., Horvath, A., Horwath, M., Khan, S., Kjeldsen, K. K., Konrad, H., Langen, P. L., Lecavalier, B., Loomis, B., Luthcke, S., McMillan, M., Melini, D., Mernild, S., Mohajerani, Y., Moore, P., Mottram, R., Mouginot, J., Moyano, G., Muir, A., Nagler, T., Nield, G., Nilsson, J., Noël, B., Otosaka, I., Pattle, M. E., Peltier, W. R., Pie, N., Rietbroek, R., Rott, H., Sørensen, L. S., Sasgen, I., Save, H., Scheuchl, B., Schrama, E., Schröder, L., Seo, K.-W., Simonsen, S. B., Slater, T., Spada, G., Sutterley, T., Talpe, M., Tarasov, L., Jan van de Berg, W., van der Wal, W., van Wessem, M., Vishwakarma, B. D., Wiese, D., Wilton, D., Wagner, T., Wouters, B., Wuite, J., and Team, T. I.: Mass balance of the Greenland Ice Sheet from 1992 to 2018, Nature, 579, 233–239,, 2019.  a, b

Slater, D. A., Straneo, F., Felikson, D., Little, C. M., Goelzer, H., Fettweis, X., and Holte, J.: Estimating Greenland tidewater glacier retreat driven by submarine melting, The Cryosphere, 13, 2489–2509,, 2019. a, b

Slater, D. A., Felikson, D., Straneo, F., Goelzer, H., Little, C. M., Morlighem, M., Fettweis, X., and Nowicki, S.: Twenty-first century ocean forcing of the Greenland ice sheet for modelling of sea level contribution , The Cryosphere, 14, 985–1008,, 2020. a, b

Sommers, A., Rajaram, H., and Morlighem, M.: SHAKTI: Subglacial Hydrology and Kinetic, Transient Interactions v1.0, Geosci. Model Dev., 11, 2955–2974,, 2018. a

Stearns, L. A. and van der Veen, C. J.: Friction at the bed does not control fast glacier flow, Science, 361, 273–277,, 2018. a

Weertman, J.: On the Sliding of Glaciers, J. Glaciol., 3, 33–38,, 1957. a

Werder, M. A., Hewitt, I. J., Schoof, C. G., and Flowers, G. E.: Modeling channelized and distributed subglacial drainage in two dimensions, J. Geophys. Res.-Earth, 118, 2140–2158,, 2013. a

Wilson, N., Straneo, F., and Heimbach, P.: Satellite-derived submarine melt rates and mass balance (2011–2015) for Greenland's largest remaining ice tongues, The Cryosphere, 11, 2773–2782,, 2017. a

Short summary
Estimates of future sea-level contribution from the Greenland ice sheet have a large uncertainty based on different origins. We conduct numerical experiments to test the sensitivity of Greenland ice sheet projections to spatial resolution. Simulations with a higher resolution unveil up to 5 % more sea-level rise compared to coarser resolutions. The sensitivity depends on the magnitude of outlet glacier retreat. When no retreat is enforced, the sensitivity exhibits an inverse behaviour.