Sensitivity of Greenland ice sheet projections to spatial resolution in higher-order simulations: the AWI contribution to ISMIP6-Greenland using ISSM

. 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 ﬂow 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 on centennial sea-level contributions from the Greenland ice sheet and study the mechanism at play. To this end, we employ the ﬁnite-element higher- 5 order ice ﬂow 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 GCM MIROC5 under the high emission scenario RCP8.5 and consider both atmospheric and oceanic forcing in full and separate scenarios. Under the full scenarios, ﬁner simulations unveil up to ∼ 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 10 forcing, the resolution dependency exhibits an opposite behaviour with about ∼ 5% more sea-level contribution in the coarser resolution. The sea-level contribution indicates a converging behaviour ≤ 1km 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 15 speeds under scenarios with outlet glacier retreat.

. Results from the initMIP-Greenland exercise (Goelzer et al., 2018). Sea-level contribution versus (   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 .

Overview of experiments 140
The ISMIP6 protocol requests the initialization mode prior or to 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:

145
initialization: experiment to retrieve the initial state of the model.
ctrl: experiment where the climate is held constant to the reference climate (from January 1991 to end of December 2100).
-ctrl_proj: experiment where the climate is held constant to the reference climate (from January 2015 to end of December 2100).

150
historical: experiment to bring model from the initialization state to ISMIP6 projection start date (from January 1991 to end of December 2014).
projection: future climate scenario (from January 2015 to end of December 2100).

Initialization experiment
The initialization state of ISSM is based on data assimilation and inversion for determining the basal friction coefficient. Before 155 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 one 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 paleoclimatic thermal spin-up and 160 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 3D temperature and watercontent 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) 165 and Goelzer et al. (2018, see submissions AWI-ISSM1 and 2) that the temperature field and its change has a negligible effect on century time-scale projections of the GrIS.
The main ingredient to the initialization is the inversion to infer the basal friction coefficient k 2 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 170 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 MEaSURE project (Joughin et al., 2016(Joughin et al., , 2018, as the data set with almost no gaps over 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 MEaSURE velocity data set (temporal coverage from 2014 to 2018). However, we ignore the contempo-175 raneity requirement in the inversion approach and put more weight on to start 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).

Historical and control experiments
In both control experiments (ctrl and ctrl_proj), the SMB and ice sheet mask remains unchanged to the reference year according 180 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  for the years from 1990 to 2015. For simplicity, the ice sheet extent remains unchanged to 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.

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 , Nowicki et al. (2020b), and Slater et al. (2019. As we aim to study the effect of grid resolution on ice mass changes, we run the future scenarios based on climate data 190 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 as 'full' in the following (RCP8.5-Rlow/med/high). In addition, we perform simulations where either the atmospheric forcing (RCP8.5-Rnone) or the marine-terminating outlet glacier retreat (OO-Rmed/high) is at play. The conducted projection 195 experiments and the corresponding experiment labels used in this study are summarized in Tab. 2 and are explained in the following sections.

Atmospheric forcing
ISMIP6 provide 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 200 to capture 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 with the anomaly defined as where SMB(x, y, t) GCM−MAR is the direct output of MAR using the GCM climate data and SMB(x, y) 1960−1989 GCM−MAR the corresponding mean value over the reference period (from January 1960 to December 1989). As the reference SMB field SMB ref (x, y), we choose the downscaled RACMO2.3p2 product  whereby a model output was averaged for the period . This period is chosen as the ice sheet is assumed close to 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 210 processes due to the evolving surface in the projection experiments are not captured. The SMB height-elevation feedback is considered with a dynamic correction SMB dyn to the SMB clim following Franco et al. (2012) SMB dyn (x, y, t) = dSMBdz(x, y, t) × (z s (x, y, t) − z ref (x, y)). (4) The surface elevation changes are taken from the ISM elevation z s (x, y, t) while running the simulation and the corresponding ISM reference elevation z ref (x, y) from the initialization state. The yearly patterns of ∆SMB(x, y, t) and dSMBdz(x, y, t) are 215 provided by ISMIP6. A cumulative SMB anomaly over the projection period is shown in Fig. 4a.

Oceanic forcing
For the oceanic forcing we rely on the empirically derived outlet glacier parametrization retreat by Slater et al. (2019Slater et al. ( , 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 the calving front, retreat and advance of marine-terminating outlet glaciers is 220 directly prescribed as a yearly series of ice front positions. (i.e., is not a result of ice velocity at calving front, calving rate and frontal melt that is used to simulate the calving front position). Here, the binary retreat masks (i.e., ice 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 submarine melting taking into account changes in ocean temperature and surface meltwater runoff from a GCM. The parametrization is not applied to the individual glaciers but 225 over a predefined geographical region. Based on the numerous retreat trajectories, a medium retreat scenario as the trajectory with the median retreat at 2100 is defined. To cover uncertainty by this approach, a low and high retreat scenarios is defined as the trajectories with the 25th and 75th percentile retreats at 2100. In the following, these retreat scenarios are termed Rlow, Rmed and Rhigh (Tab. 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 experience no retreat of marine-terminating outlet glaciers.

Comparability of experiments
A central question about resolution dependence is always "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 to observational targets or to achieve a reduced model drift. However, we decided here to keep model parameters (e.g. inversion 235 parameters) and parameterizations (e.g. sub-grid scheme at grounding line) unchanged for all simulations. Similarly for the retreat masks, we rely on binary retreat masks for all adopted resolutions although the ISMIP6 protocol requests a sub-grid scheme for coarse resolution models. On the 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 3 Results

Historical scenario
To evaluate the modelling decisions pertaining 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 MEaSURE data sets (250 m grid spacing) for ice thickness and surface velocity, respectively. As these data are used in the data assimilation and inversion, velocity and thickness are not independent 250 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 at the ISMIP6 projection start date.
Since the results are qualitatively similar for each grid simulation (Figs. S1, S2 and S3), 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 The stored volumes, ice extent and spatially integrated SMB is among all grid simulations rather similar (V = 7.28 m SLE± 0.2%, A = 1.787 × 10 6 km 2 ± 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

275
intrinsically simulated ice discharge Q yields 207 Gt a −1 to 341 Gt a −1 for the G4000 and G750 simulations, respectively.

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×10 8 km 2 . In the following, the In the ctrl experiment, the transient response (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 with a low rate of volume change. As the initialization states are presumably different across the employed grids, we expect a 285 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 290  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 305 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 mm 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 inline with earlier studies where this effect was already 310 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/med scenarios exhibits a spread of 10.7 mm and 13.6 mm, respectively, which is in the order of the absolute magnitude. A notable feature for all conducted simulations is, that 315 the sea-level contribution in each individual experiment converges with increasing resolution. .5-Rhigh show a peak in mass loss at the finest resolution, whereby a peak in mass loss is attained in the G2000 simulation for RCP8.5-Rlow. For the latter, it is worth to mention 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 get a more in-depth insight into the role of atmospheric and oceanic forcing on grid-resolution. It is worth to mention that the qualitative behaviour of the detected grid-dependent mass loss remains similar when the projections are corrected with the ctrl_proj experiment (Tab. 3).

Mass partitioning 330
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 In the full experiments RCP8.5-Rlow/med/high, an increase in resolution enhances the importance of dynamic contribution.

345
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 350 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 explains 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 355 finer resolutions. This finding corroborates with the grid-dependent increase of 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.

Outlet glacier response 360
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 the retreat, or other non obvious factors. We conduct an in-depth analysis of numerous prominent outlet glaciers at GrIS (Fig. S6 and table with analysis provided as separate SI). The responses of most of the outlet glacier reveal the deduced grid-dependent behaviour where higher resolutions cause an enhanced discharge. This is exemplary illustrated in Fig. 11a for Helheim Glacier. However, this behaviour could not be 365 adopted to 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 from the finer resolution. This topographic effect is restricted to narrow confined outlet glaciers that obey a characteristic width in 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. Theses glaciers seem to have a qualitatively equal behaviour for glacier speed-up 370 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.
Glaciers that are converted from a marine-terminating to a land-terminating glacier by retreating out of the water build an own class. These glaciers are no longer subject of the retreat and show a collapse in ice discharge regardless of the grid 375 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 employed retreat parametrization that lacks information of the bedrock topography, such as topographic highs and lows.

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

Discussion
The simulations presented here show that the projected sea-level contribution is sensitive to the spatial resolution. The sensitiv-385 ity effect depends on the climate forcing, with oceanic and atmospheric forcings showing opposite and non-trivial responses.
The simulations have turned out 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, dependent 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 390 explains the different responses. 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.
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 effec-395 tive pressure (Eq. 1), basal sliding velocities increase stronger in the finer resolution. This feedback is enhanced as the SMB perturbations lead to a decrease in ice thickness, hence in a decrease of 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 time scales as the surface melt is amplified by the positive surface mass balance-elevation feedback exposing the ice surface to higher air temperatures.

400
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 slow-down. In the corresponding G750 simulation most of the outlet glaciers show a speed-up but this effect is very localized and do not reach far upstream. Therefore, we conclude, that the pronounced 405 decrease of the effective pressure along with the acceleration of outlet glacier is a dominant mechanism controlling the griddependent spread.
In order to investigate whether the response behaviour is an effect by purely reducing the grid size, we repeated the OO-Rhigh and RCP8.5-Rhigh experiments with a G1000 simulation using re-gridded bedrock topography and friction coefficient from the G4000 initial state (simulations are not shown). This setup adopts a high-resolution grid but omits detailed information from the 410 high-resolution input data. Projected mass loss by this setups is closer to the G4000 simulation. They, therefore, demonstrate that a high model resolution alone is insufficient to explain the grid-dependent sea-level contribution. As a consequence, a driving cause for 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.
A convergence of the grid-dependent estimates of sea-level contribution emerges around RES high ≤ 1 km. This value corroborates with 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 whether a model res-440 olution 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 parametrization is insensitive to bed undulations, resolving the outlet glacier cross-sections is important for accurately model ice discharge. Since the glacier cross-sections are reasonable well approximated in G750, we do not expect that a resolving the geometry higher would drastically alter ice discharge rates. On the other hand, there are indications that at a resolution of 750 m the HO solution 445 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 seem 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. That might indicate that the underlying processes are not fully converged and still causing changes in mass loss trends.

450
Our grid-dependent results under atmospheric only forcing correlates with the finding in Goelzer et al. (2018, Fig . 1) and the Exp. C2 in Greve and Herzfeld (2013, Figs. 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., 2019), whereby in our study the effect is attributed to the change is ice discharge of the ice sheet. The cause for the griddependent behaviour in Greve and Herzfeld (2013) is not specified further. Still, it is worth mentioning that they report a much 455 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 scenarios 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).

460
Besides the fixed calving front in the atmospheric only scenario, further limitations of our study must be noted.  Leguy et al., 2014;Joughin et al., 2019). In brief, these type 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 time scales.
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 re-470 sults in an overall comparable pattern for the flow velocities ( Fig. S1 and S2) and basal friction coefficient (Fig. S3) for all grids. However, on smaller scales, the inversion approach produces significantly different k 2 in many glacier basins. Recalling the relationship between N and k 2 , these different patterns are plausible but could potentially be a result from non optimal inversion parameters. However, this different spatial patterns might 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 475 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 comparable magnitude as 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% . A feedback that we have not considered is the enhanced surface melt influencing the basal conditions. Given that we greatly simplify the 480 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. Subglacial hydrology model has shown the localized effect on N , which is likely having consequences on the spread between the employed grid resolutions (e.g. Werder et al., 2013;de Fleurian et al., 2016;Sommers et al., 2018;Beyer et al., 2018).

485
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 increases the widespread inland flow acceleration and enhances the rate of mass loss. However, the change in τ b and τ d is very pronounced in and around the main trunks and quite differently for the 490 adopted grids (see Fig. S20 and 21). These patterns exemplify the need for resolving the shear margins particularly high, which we have not fully accomplished in this study as shear margins are becoming, in numerous cases, sub-grid phenomena. This may be a reason for under-representing glacier velocities inside the main trunk and over-estimate 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).

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 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-500 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 recognise for ice sheet models that have SMB as the dominant mass loss driver.
bedrock undulations experience a similar response in all model resolution. 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 glacier. Since we rely on a greatly simplified effective pressure parametrization, further work with 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 grids and physics adequately.