The effect of overshooting 1.5 °C global warming on the mass loss of the Greenland ice sheet

Sea-level rise associated with changing climate is expected to pose a major challenge for societies. Based on the efforts of COP21 to limit global warming to 2.0 ^∘C or even 1.5 ^∘C by the end of the 21st century (Paris Agreement), we simulate the future contribution of the Greenland ice sheet (GrIS) to sea-level change under the low emission Representative Concentration Pathway (RCP) 2.6 scenario. The Ice Sheet System Model (ISSM) with higher-order approximation is used and initialized with a hybrid approach of spin-up and data assimilation. For three general circulation models (GCMs: HadGEM2-ES, IPSL-CM5A-LR, MIROC5) the projections are conducted up to 2300 with forcing fields for surface mass balance (SMB) and ice surface temperature (T_s) computed by the surface energy balance model of intermediate complexity (SEMIC). The projected sea-level rise ranges between 21–38 mm by 2100 and 36–85 mm by 2300. According to the three GCMs used, global warming will exceed 1.5 ^∘C early in the 21st century. The RCP2.6 peak and decline scenario is therefore manually adjusted in another set of experiments to suppress the 1.5 ^∘C overshooting effect. These scenarios show a sea-level contribution that is on average about 38 % and 31 % less by 2100 and 2300, respectively. For some experiments, the rate of mass loss in the 23rd century does not exclude a stable ice sheet in the future. This is due to a spatially integrated SMB that remains positive and reaches values similar to the present day in the latter half of the simulation period. Although the mean SMB is reduced in the warmer climate, a future steady-state ice sheet with lower surface elevation and hence volume might be possible. Our results indicate that uncertainties in the projections stem from the underlying GCM climate data used to calculate the surface mass balance. However, the RCP2.6 scenario will lead to significant changes in the GrIS, including elevation changes of up to 100 m. The sea-level contribution estimated in this study may serve as a lower bound for the RCP2.6 scenario, as the currently observed sea-level rise is not reached in any of the experiments; this is attributed to processes (e.g. ocean forcing) not yet represented by the model, but proven to play a major role in GrIS mass loss.

1 Introduction

Within the past decade, the Greenland ice sheet has contributed about 20 % to sea-level rise (Rietbroek et al., 2016), caused by the acceleration of outlet glaciers and changes in the surface mass balance (Enderlin et al., 2014). In the past decades, these changes in surface mass balance contributed about 60 % of the ice sheet mass loss, whereas 40 % is attributed to increasing discharge (van den Broeke et al., 2016). The question arises of what impact the GrIS will have on sea-level change in the next decades and centuries.

Negotiated during COP21, the Paris Agreement's aim is to keep the global temperature rise in this century well below 2 ^∘C above pre-industrial levels and to pursue efforts to limit the temperature increase even further to 1.5  ^∘C (UNFCCC, 2015). However, the statement of holding global temperature below 2 ^∘C implies keeping global warming below the 2 ^∘C limit over the course of the entire century and beyond, while efforts to limit the temperature increase to 1.5 ^∘C are often interpreted as allowing for a potential overshoot before returning to below 1.5 ^∘C (Rogelj et al., 2015). Here we selected the Representative Concentration Pathway (Moss et al., 2010) 2.6 as the lowest emission scenario considered within CMIP5 and in line with a 1.5 or 2 ^∘C limit on global warming. Depending on the general circulation models (GCMs) considered, the global temperature change over time varies considerably, although the political target is met by 2100. Whereas some models in RCP2.6 do not exceed the limit of 1.5 or 2.0 ^∘C global warming before 2100, other models do and exhibit subsequent cooling (Frieler et al., 2017).

While global temperature rise may be limited to 1.5 or 2 ^∘C by 2100, warming over Greenland is enhanced due to the Arctic amplification (Pithan and Mauritsen, 2014), has already exceeded 1.5 ^∘C (relative to 1951–1980) in the past decade (GISTEMP Team, 2018), and may exceed 4 ^∘C by 2100. This yields more than 2 ^∘C warming by 2100 and could therefore have a considerable impact on ice sheet mass loss over Greenland. This implies an enlargement of the ablation zone and goes along with a decline in SMB. However, it is currently unclear how fast the GrIS could react to cooling and recovery of SMB, as ice sheets also react dynamically to atmospheric forcing.

Recent large-scale ice sheet modelling attempts to project the contribution of the GrIS under RCP2.6 warming scenarios are very scarce. Fürst et al. (2015) conducted an extensive study to simulate future ice volume changes driven by both atmospheric and oceanic temperature changes for all four RCP scenarios. For the RCP2.6 scenario they estimate a sea-level contribution of 42.3±18.0 mm by 2100 and 88.2±44.8 mm by 2300. The value by 2100 is in line with estimates given by the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC AR5, IPCC, 2013). The AR5 range for RCP2.6 lies between 10 and 100 mm by 2100 (the value depends on whether ice dynamical feedbacks are considered or not).

The GrIS response to projections of future climate change are usually studied with a numerical ice sheet model (ISM) forced with climate data. ISM response is subject to the dynamical part and the surface mass balance (SMB). In the past, ISMs often used the rather simple and empirically based positive degree day (PDD) scheme, in which the PDD index is used to compute melt, run-off, and ice surface temperature from atmospheric temperature and precipitation (Huybrechts et al., 1991). One disadvantage of the PDD method is that the involved PDD parameters are tuned to correctly represent present-day melting rates but may fail to represent past or future climates (Bauer and Ganopolski, 2017; Bougamont et al., 2007). On one far end of model complexity, a regional climate model (RCM) resolves most processes at the ice–atmosphere interface and in the upper firn layers, such as RACMO (Noël et al., 2018) or MAR (Fettweis et al., 2017), with higher spatial and temporal resolution than GCMs. RCMs have been shown to be quite successful in reproducing the current SMB of the GrIS. However, as they are computationally expensive, an intermediate way would be more efficient to balance the computational costs and parameterization of processes, such as the surface energy balance model of intermediate complexity (Krapp et al., 2017), which is employed in this study.

Here we target RCP2.6 peak and decline scenarios in particular to study the GrIS response to overshooting with a numerical ISM. The projections are driven with climate data output from the CMIP5 RCP2.6 scenario provided by the ISIMIP2b project for different GCMs (Frieler et al., 2017). To obtain ice surface temperature and surface mass balance from the atmospheric fields, the surface energy balance model SEMIC (Krapp et al., 2017) is applied. The SEMIC (Sect. 2.1) is driven off-line to the ISM and therefore the climate forcing is one-way coupled and applied as anomalies to the ISM. The advantage of this one-way coupling is the lower computational cost, allowing for reasonably high spatial and temporal resolution of the ISM. In order to study the effect of overshooting, we design an RCP2.6-like scenario without overshoot by manually stabilizing the forcing at 1.5 ^∘C.

For modelling the flow dynamics and future evolution of the GrIS under RCP2.6 scenarios, the thermo-mechanical coupled Ice Sheet System Model (Larour et al., 2012) with a Blatter–Pattyn-type higher-order momentum balance (Blatter, 1995; Pattyn, 2003) is applied (Sect. 2.5). A crucial prerequisite for projections is a reasonable initial state of the ice sheet in terms of ice thickness, ice extent, and ice velocities. Beside starting projections with the most realistic setting, the prevention of a model shock after switching from the initialization procedure to projections is very important. Both have been a major issue in the past, which gave rise to an international benchmark experiment, initMIP Greenland (Goelzer et al., 2018), for finding optimal strategies to derive initial states for the ice velocity and temperature fields. Here, we apply a hybrid approach between a thermal paleo-spin-up and data assimilation.

Before driving the projections, the SMB forcing is validated thoroughly against RACMO. Then we explore the response of the GrIS and its contribution to sea-level rise under the RCP2.6 scenario and a modified RCP2.6 scenario without overshoot.

2 Model description

2.1 Energy balance model

Thermo-mechanical ISMs require the annual mean surface temperatures and annual mean surface mass balance of ice as boundary conditions at the surface. To derive these ice-sheet-specific quantities, we use the surface energy balance model of intermediate complexity (Krapp et al., 2017). Although we only apply SEMIC and do not adjust the parameters of SEMIC, SEMIC is described briefly. SEMIC computes the mass and energy balance of snow and/or ice surface. In order to tune parameters for a number of processes, Krapp et al. (2017) performed an optimization for the GrIS forced with regional climate model data (MAR). These parameters have been used in our study, too. The energy balance equation reads as

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }{c}_{\mathrm{eff}}{\displaystyle \frac{\mathrm{d}{T}_{\mathrm{s}}}{\mathrm{d}t}}=(\mathrm{1}-\mathit{\alpha })\cdot {\mathrm{SW}}^{\downarrow }-{\mathrm{LW}}^{\uparrow }+{\mathrm{LW}}^{\downarrow }-H_{\mathrm{S}}-H_{\mathrm{L}}-Q_{\mathrm{M}/\mathrm{R}},\end{array}$$

where α is the surface albedo that is parameterized with the snow height (Oerlemans and Knap, 1998). The downwelling shortwave SW^↓ and downwelling longwave radiation LW^↓ at the surface are provided as atmospheric forcing (Sect. 2.2). The upwelling longwave radiation LW^↑ is described by the Stefan–Boltzmann law. The latent H_L and sensible H_S heat fluxes are estimated by the respective bulk approach (Gill, 1982). The residual heat flux Q_M∕R is calculated from the difference of melting M and refreezing R and keeps track of any heat flux surplus or deficit to keep the ice surface temperature T_s below or equal to 0 ^∘C over snow and ice.

The surface mass balance (SMB) in SEMIC is considered as follows:

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }\mathrm{SMB}=P_{\mathrm{s}}-\mathrm{SU}-M-R,\end{array}$$

where P_s is the rate of snowfall and SU the sublimation rate, which is directly related to the latent heat flux. The melt rate M is dependent on the snow height; if all snow has melted down, the excess energy is used to melt the underlying ice. The refreezing R is calculated differently for available meltwater or rainfall. Moreover, the porous snowpack could retain a limited amount of meltwater, while over ice surfaces refreezing is neglected and all melted ice is treated as run-off. In SEMIC, the total melt rate and refreezing rate are calculated from the available energy during the course of 1 day. As the set of equations is solved using an explicit time step scheme with a time step of 1 day, a parameterization for the diurnal cycle (a cosine function) accounts for thawing and freezing over a day. This reduces complexity; the one-layer snowpack model saves computational time and allows for integration on multi-millennial timescales as opposed to more sophisticated multilayer snowpack models. Further details are given by Krapp et al. (2017).

2.2 Atmospheric forcing

Here we targeted peak and decline scenarios in particular, temporarily exceeding a given temperature limit of global warming to 2.0 ^∘C or even 1.5 ^∘C by the end of 2100. From the official extended RCP2.6 scenarios (Meinshausen et al., 2011), we have selected GCMs which cover the CMIP5 historical scenario, the RCP2.6 scenario until 2299, and reveal an overshoot in annual global mean near-surface temperature change relative to pre-industrial levels (1661–1860). Three different GCMs were used in our study: IPSL-CM5A-LR (L'Institut Pierre-Simon Laplace Coupled Model, version 5, low resolution), MIROC5 (Model for Interdisciplinary Research on Climate, version 5), and HadGEM2-ES (Hadley Centre Global Environmental Model 2, Earth System). Instead of the full acronyms we use IPSL, HadGEM2, and MIROC5 in the following text. The GCM output was provided and prepared by the ISIMIP2b project following a strict simulation protocol (Frieler et al., 2017). Figure 1a displays the temporal evolution of the annual global mean near-surface air temperature T_a for those GCMs for the historical simulation up to 2005 continued with the RCP2.6 simulation up to 2299. Global mean temperature projections from IPSL, HadGEM2, and MIROC5 under RCP2.6 exceed 1.5 ^∘C relative to pre-industrial levels in the second half of the 21st century. While global mean temperature change returns to 1.5 ^∘C or even slightly lower by 2299 in HadGEM2, it only reaches about 2 ^∘C in IPSL by 2299. For MIROC5, temperature stabilizes at about 1.5 ^∘C during the second half of the 21st century. In order to determine the onset of overshoot we scan the historical and RCP2.6 scenarios of the individual GCMs to identify the time when the global warming reaches 1.5 ^∘C in a 30-year moving window above pre-industrial levels. The characteristic date of overshooting 1.5 ^∘C for HadGEM2 is by 2023; MIROC5 reaches this level by 2043, while IPSL reaches this point by 2009 (coloured dots in Fig. 1).

Figure 1 Time series of annual global mean near-surface temperature change (a) and over the GrIS (b) for all three GCMs relative to 1661–1880. The thick line is a 30-year moving mean. The coloured dots represent the onset years of overshooting 1.5 ^∘C in the global mean near-surface air temperature in a 30-year moving window relative to pre-industrial levels. The light gray shaded area indicates the reused time period for the scenario without overshoot.

The phenomenon that tends to produce a larger change in temperature near the poles was termed polar amplification. Particularly, it enhances the increase in global mean air temperature over Arctic areas (referred to here as Arctic amplification). Generally, the CMIP5 models show an annual average warming factor over the Arctic between 2.2 and 2.4 times the global average warming (IPCC, 2013). As mechanisms creating the Arctic amplification may be represented to different extents in the GCMs, the level of future amplification is different across the GrIS. The three GCMs used in this study represent this trend to differing extents over the GrIS¹ (Figs. 1 and 2). For HadGEM2 and IPSL the Arctic compared to global warming is amplified relatively similarly (warming approximately 4 ^∘C relative to 1661–1860). In contrast, MIROC5 reveals a considerably lower Arctic amplification (warming approximately 3 ^∘C relative to 1661–1860). In terms of global and Arctic future annual mean near-surface temperatures, MIROC5 offers the lowest and IPSL the highest forcing.

The ISIMIP2b atmospheric forcing data are CMIP5 climate model output data that have been spatially interpolated to a regular $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ latitude–longitude grid and bias-corrected using the observational dataset EWEMBI (Frieler et al., 2017; Lange, 2018). To drive the SEMIC, we need to provide the atmospheric forcing (consisting of incoming shortwave radiation SW^↓, longwave radiation LW^↓, near-surface air temperature T_a, surface wind speed u_s, near-surface specific humidity q_a, surface air pressure p_s, snowfall rate P_s, and rainfall rate P_r). These fields are available from the output data of the three GCMs. SEMIC is driven by the daily input of the GCMs, while the output is the cumulative surface mass balance and the mean surface temperature over each year.

Erokhina et al. (2017)Vizcaíno et al. (2010)Vizcaíno et al. (2010)Vizcaíno et al. (2010)

Table 1 Lapse rates and height–desertification relationship for initial corrections of the GCM output fields near-surface air temperature T_a, precipitation of snow P_s, precipitation of rain P_r, and downwelling longwave radiation LW^↓ used as input for SEMIC. Here, h^ref=2000 m is the reference height and ${\mathit{\gamma }}_p=-\mathrm{0.6931}$ km⁻¹ is the desertification coefficient.

Figure 2 Scatter plot of annual mean near-surface air temperature change relative to pre-industrial levels over GrIS versus annual global mean near-surface air temperature change for the years 1861–2299. The gray line is the identity line.

Given the differences in resolution between the GCMs and ISSM, a vertical downscaling procedure is applied to the atmospheric forcing fields. First, the atmospheric fields are interpolated bilinearly from the GCM grid onto a regular high-resolution 0.05^∘ grid on which SEMIC is run. As a result, the output fields of SEMIC are conservatively interpolated onto the unstructured ISSM grid. This two-step procedure is not necessary, but currently the easiest way from a technical standpoint. For future applications we will avoid the intermediate interpolation and run SEMIC directly on the unstructured target ISSM grid. To account for the difference in ice sheet surface topography between GCMs and ISSM, corrections for several quantities (⋅) denoted by (⋅)^cor are initially performed. We follow the corrections proposed by Vizcaíno et al. (2010),

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }(\cdot {)}^{\mathrm{cor}}=\left(h_{\mathrm{s}}^{\mathrm{SEMIC}}-h_{\mathrm{s}}^{\mathrm{GCM}}\right){\mathit{\gamma }}_{(\cdot )},\end{array}$$

with the lapse rates γ_(⋅) shown in Table 1 and $h_{\mathrm{s}}^{\mathrm{SEMIC}}$ equal to the ISSM ice surface elevation at the initial state. Subsequently, SEMIC computes the ice surface temperature T_s and the SMB based on these corrected input values.

SEMIC is applied as developed by Krapp et al. (2017). These authors perform a particle-swarm optimization to calibrate model parameters for the GrIS and validate them against the RCM MAR. We adopt their derived parameters here. The parameter tuning aimed to find a parameter set which gives the best fit between SMB and ice temperature T_s of SEMIC with only a limited number of processes and simpler parameterizations compared to a more complex RCM. An RCM is typically validated against reanalysis data and observations; therefore, we assume the tuned parameters are most reliable to represent the processes and parameterizations within SEMIC. In terms of process description, the optimized SEMIC configuration leads to the best possible SMB and T_s fields when MAR is used as forcing. If SEMIC is tuned with another RCM (e.g. RACMO or HIRHAM), the parameters will be different. Here, a separate tuning for each GCM would be required due to the differences (e.g. the timing of maximum warming, the length of an overshoot) among the GCMs used in this study. This basically means compensating for too-low near-surface temperatures, for example, with SEMIC parameters, which would offset the whole comparison of GCM forcing. Furthermore, these additional tuning steps would make the benefit of having a semi-complexity model with low costs meaningless.

Figure 3 compares averaged SMB fields for the time period 1960–1990 from RACMO2.3 and, as an example, the SMB derived by forcing SEMIC with HadGEM2. The pattern of the SMB derived by forcing SEMIC with IPSL or MIROC5 is the same as when using HadGEM2. The comparison in Fig. 3 shows that the large-scale patterns of the SMB fields are in fairly good agreement, while the small-scale pattern and magnitudes of the GCM-based SMB do not match the RACMO2.3 SMB. Although the coarse GCM-based forcing has undergone a downscaling of particular fields and is processed in SEMIC with a higher resolution, the atmospheric fields over the ice sheet still lack details compared to an RCM. This is due to the fact that the forcing from a GCM implies different characteristics, like smoother gradients and a less-resolved geometry, compared to the RCM. The direct output of the SMB from SEMIC to the RACMO2.3 field has a misfit of about 2 m a⁻¹ and a correlation coefficient of R²=0.5. Additionally, the spatially integrated SMB for the averaged time period differs by up to 200 Gt a⁻¹. For HadGEM2, IPSL, and MIROC5 the values are 536, 496, and 614 Gt a⁻¹, respectively. In contrast, the corresponding value for RACMO2.3 is 403 Gt a⁻¹. Therefore, we conclude that the absolute fields from SEMIC are not ideal for our purpose.

Figure 3 Comparison of surface mass balance fields averaged for the time period 1960–1990; (a) surface mass balance derived by forcing SEMIC with climate data from HadGEM2; (b) surface mass balance of RACMO2.3 (Noël et al., 2018).

Instead of using the absolute SEMIC output fields (SMB and T_s) directly to force the numerical ice flow model ISSM, we rely on an anomaly method. The climatic boundary conditions applied here consist of a reference field onto which climate change anomalies from SEMIC are superimposed. The initialization of the ice flow model based on data assimilation (Sect. 2.6 below) makes it possible to use forcing data from high-resolution RCMs that were run on the same ice sheet mask and ice surface topography. As the reference SMB field we choose the downscaled RACMO2.3 product (Noël et al., 2018) whereby a model output was averaged for the time period 1960–1990, denoted $\stackrel{\mathrm{‾}}{\mathrm{SMB}}$(1960–1990)_RACMO. The reference period 1960–1990 is chosen as the ice sheet is assumed close to steady state in this period (Ettema et al., 2009). The climatic SMB that is used as future climate forcing reads

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }{\mathrm{SMB}}_{\mathrm{clim}}(x,y,t)={\stackrel{\mathrm{‾}}{\mathrm{SMB}}}_{\mathrm{RACMO}}^{(\mathrm{1960}-\mathrm{1990})}(x,y)+\mathrm{\Delta }\mathrm{SMB}(x,y,t),\end{array}$$

with the anomaly defined as

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }\mathrm{\Delta }\mathrm{SMB}(x,y,t)={\mathrm{SMB}}_{\mathrm{SEMIC}}(x,y,t)-{\stackrel{\mathrm{‾}}{\mathrm{SMB}}}_{\mathrm{SEMIC}}^{(\mathrm{1960}-\mathrm{1990})}(x,y),\end{array}$$

where $t=\mathit{\{}\mathrm{1960},\mathrm{1961},\mathrm{\dots },\mathrm{2299}\mathit{\}}$. Note that the historical scenario is run from 1960–2005 and followed by the RCP2.6 scenario from 2006–2299. In an ideal case, both reference terms $\stackrel{\mathrm{‾}}{\mathrm{SMB}}$(1960–1990)_RACMO and $\stackrel{\mathrm{‾}}{\mathrm{SMB}}$(1960–1990)_SEMIC will cancel out and the absolute climatic forcing ${\mathrm{SMB}}_{\mathrm{SEMIC}}(x,y,t)$ would remain. This is certainly not the case and the equation must be interpreted as having the RACMO2.3 reference field (with a good spatial distribution) as a background field with the trends from SEMIC superimposed.

The same equations hold for the temperature imposed on the ice surface. This ensures that the unforced control experiment produces identical behaviour for each GCM. Results for future projections depend only on the atmospheric GCM input or, similarly, SEMIC output and therefore the results can be compared quantitatively. In the following text, the constructed SMB fields according to Eq. (4) are referred to as SEMIC-HadGEM2, SEMIC-IPSL, SEMIC-MIROC5, or in general as SEMIC-GCM.

In the presented study, the ice flow model is forced with the off-line-processed SEMIC output. As the ice sheet evolves in response to climate change, local climate feedback processes are not captured. Most importantly, this includes the interaction of the ice surface between air temperature and precipitation, which in turn affects the surface mass balance. The SMB feedback process is considered with a dynamic correction to the SMB_clim (see Sect. 2.7 below). This correction is applied within ISSM and to the surface mass balance term only.

Figure 4 Time series of the annual mean integrated SMB_clim (Gt yr⁻¹) according to Eq. (4) for all three SEMIC-GCMs under RCP2.6 forcing (a) and RCP2.6 forcing without overshoot (b). The solid line is a 30-year and 15-year moving mean in (a) and (b), respectively. The darker gray shading and black line represent the range and mean of SMB between 1981 and 2010 from Polar Portal (http://polarportal.dk/forsiden/, last access: 8 October 2018). The dashed line shows the SMB time series of RACMO2.3 (Noël et al., 2018) from 1958–2016. The coloured dots represent the onset years of overshooting 1.5 ^∘C in the global mean near-surface air temperature in a 30-year moving window relative to pre-industrial levels. The light gray shaded time period indicates the repeated SMB forcing taken from the RCP2.6 scenario for the scenario without overshoot.

2.3 Modified RCP2.6 scenario without overshoot

The global climate warming of the selected GCMs exceeds the political target of 1.5 ^∘C during the 20th century, although RCP2.6 is the strongest mitigation scenario (Moss et al., 2010). In order to estimate the effect of overshooting on the projected sea-level contribution from the GrIS we manually construct an RCP2.6-like scenario without an overshoot, assuming an immediate climate stabilization at that time when 1.5 ^∘C is reached. The characteristic time of overshooting 1.5 ^∘C for HadGEM2 is by 2023; MIROC5 reaches this level by 2043, while IPSL reaches this point by 2009. Before reaching this threshold, the unaltered historical and RCP2.6 forcing is applied. The extension of the forcing from these characteristic times is of crucial importance. To avoid an unphysical step change, the climate in the repeated time period should stabilize (i.e. no long-term trends in temperature change) close to 1.5 ^∘C warming. In order to account for decadal variability, i.e. extreme years, we reuse the climatic forcing fields from 2250–2280 until the end of the simulation (light gray shaded areas in Figs. 1 and 4). In this time window, the warming is close to 1.5 ^∘C and exhibits a frequent number of extreme years. Other time windows might also be feasible (e.g. the last 30 or 50 years), but will likely not change the forcing substantially. In the following, the modified RCP2.6-like scenario without overshoot is termed RCP2.6 without overshoot.

2.4 Assessment of SMB forcing

We want to emphasize the fact that we do not intend to validate the energy balance model SEMIC itself, but rather assess whether the obtained SMB fields by forcing SEMIC with the GCMs are plausible. In order to do so, the obtained climatic SMB_clim (Eq. 4), the resulting SMB patterns, and time series are compared with other available datasets. Beside the spatial pattern of the surface mass balance, the time series of the SMB over Greenland illustrates what the ice sheet's total surface gains and losses have been. The constructed SMB forcings for the RCP2.6 scenario with and without overshoot are shown in Fig. 4a and b, respectively. The gray shaded box and black line depict the range and the mean SMB between 1981 and 2010 from Polar Portal (http://polarportal.dk/forsiden/, last access: 8 October 2018) derived from a combination of observations and a weather model for Greenland (HIRLAM “newsnow”). The dashed black line shows the results from the RACMO2.3 product. The spatially integrated SMB magnitude of each SEMIC-GCM is consistent with the RACMO2.3 and Polar Portal data. The drop in SMB after 2000 is present in all three SEMIC-GCMs and RACMO2.3.

Table 2 Annual mean integrated SMB (Gt yr⁻¹) covering various periods. Time series of SMB_clim for the SEMIC-GCMs are calculated by Eq. (4) for the RCP2.6 scenario. The column “1.5 ^∘C reached” shows a 30-year mean at the characteristic time of overshooting 1.5 ^∘C. The anomaly in SMB (ΔSMB) is in 2080–2099 with respect to 1980–1999.

^* MAR forced with GCM NorESM1-M under the RCP2.6 scenario (Fettweis et al., 2013).

Figure 5 Inter-annual SMB variability for all SEMIC-GCMs (coloured lines) and RACMO2.3 (black line) calculated from consecutive years; $\mathrm{\Delta }\mathrm{SMB}={\mathrm{SMB}}_t-{\mathrm{SMB}}_{t-\mathrm{1}}$. The thick lines are a 30-year moving mean calculated from the yearly data (thin lines).

For SEMIC-HadGEM2 the spatially integrated SMB remains around 200 Gt a⁻¹ after 2050. The SMB for SEMIC-IPSL recovers from 2050 onwards and shows an increase from around 200 Gt a⁻¹ to around 350 Gt a⁻¹ by 2300. SEMIC-MIROC5 reveals the lowest SMB change over time and recovers after 2050 from 250 to 300–350 Gt a⁻¹ by 2300. By 2300, the SMB of SEMIC-IPSL and SEMIC-MIROC5 is slightly below the present-day magnitude. However, the decline of SMB for the RCP2.6 scenario roughly corresponds with MAR results forced with the GCM NorESM1-M under the RCP2.6 scenario (Fettweis et al., 2013 and last column in Table 2), although it is not strictly comparable because they use different GCM climate data. They estimated a loss of $-\mathrm{124}\pm \mathrm{100}$ Gt a⁻¹ in 2080–2099 relative to 1980–1999.

Table 2 shows the annual mean integrated SMB over the entire GrIS for various periods. Averaged over most of the periods, the annual mean integrated SMB is rather similar among the models. Most obvious are the differences between the SEMIC-GCMs for the period 1997–2016. The year 1997 was identified as the critical time of Greenland's peripheral glacier and ice cap mass balance decrease (Noël et al., 2017). For this period of declining SMB, the SEMIC-HadGEM2 agrees well with the RACMO2.3 product, while the spatially integrated SMBs for SEMIC-IPSL and SEMIC-MIROC5 are ∼40 and 50 Gt a⁻¹ larger, respectively.

For the available RACMO2.3 time series and the SEMIC-GCMs, we have computed the inter-annual SMB variability (Fig. 5). The SMB variability is similar to RACMO2.3 in terms of frequency and amplitude but is not coherent among all models because the GCMs have their own internal variability. For the time period 1960–2016, the overall surface mass balance difference over the ice sheet between SEMIC-GCM and RACMO2.3 is almost zero with −0.007, 0.016, and 0.0200 m a⁻¹ for SEMIC-HadGEM2, SEMIC-IPSL, and SEMIC-MIROC5, respectively. These numbers are in the same range as given by Krapp et al. (2017) for the comparison between SEMIC and MAR.

2.5 Ice flow model

The ice flow and thermodynamic evolution of the GrIS are approximated using the finite-element-based ISSM. The model has been applied successfully to both large ice sheets (Bindschadler et al., 2013; Goelzer et al., 2018; Nowicki et al., 2013) and is also used for studies of individual drainage basins of Greenland, e.g. the North-East Greenland Ice Stream (NEGIS) (Choi et al., 2017), Jakobshavn Isbræ (Bondzio et al., 2016, 2017), and Store Glacier (Morlighem et al., 2016). Here, we use an incompressible non-Newtonian constitutive relation with viscosity dependent on temperature, microscopic water content, and strain rate, while neglecting the softening effect of damage or impurities. The BP approximation to the Stokes momentum balance equation is employed in order to account for longitudinal and transverse stress gradients.

ISSM is specified with kinematic boundary conditions at the upper and lower boundary of the ice sheet. The upper boundary incorporates the climatic forcing obtained from SEMIC as explained above, i.e. the surface mass balance and ice surface temperature. The ice surface temperature is prescribed through Dirichlet boundary conditions. The base of grounded ice is specified as both impenetrable with the bedrock and in balance with the rate of basal melting. At the base of floating ice we use a Neumann boundary condition that parameterizes the heat flux at the ice–ocean interface (Larour et al., 2012). The basal melt rate below ice shelves is parameterized with a Beckmann–Goosse relationship (Beckmann and Goosse, 2003). The melt factor is roughly adjusted such that melting rates correspond to literature values (Wilson et al., 2017). Within this study the basal melt rate is not a focus and hence the basal melt underneath floating tongues or vertical calving fronts of tidewater glaciers are not changed. Once the pressure melting point at the grounded ice is reached, melting is calculated from basal frictional heating and the heat flux difference at the ice–bed interface. At the ice base, sliding is allowed everywhere and the basal drag, τ_b, is written using Coulomb friction:

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle }{\mathit{\tau }}_{\mathrm{b}}=-k^{\mathrm{2}}N{\mathit{v}}_{\mathrm{b}},\end{array}$$

where v_b is the basal velocity vector tangential to the glacier base and k² a constant. The effective pressure is defined as $N={\mathit{\varrho }}_{\mathrm{i}}gH+{\mathit{\varrho }}_{\mathrm{w}}g{h}_{\mathrm{b}}$, where H is the ice thickness, h_b the glacier base, and ϱ_i=910 kg m⁻³ and ϱ_w=1028 kg m⁻³ the densities for ice and seawater, respectively. We apply water pressure at the calving front of marine-terminating glaciers and observed surface velocities (Rignot and Mouginot, 2012) at the ice front of land-terminating glaciers. A traction-free boundary condition is imposed at the ice–air interface.

Geothermal heat flows into the ice in contact with bedrock and adjusts dynamically to the thermal state of the base (Aschwanden et al., 2012; Kleiner et al., 2015). The spatial pattern of the geothermal flux is taken from Greve (2005).

For all simulations, the ice front is fixed in time, and a minimum ice thickness of 10 m is applied. This implies that calving and melting exactly compensate for the outflow through the margins and initially glaciated points are not allowed to become ice free. However, regions that reach this minimum thickness have retreated. The grounding line is allowed to evolve freely according to a sub-grid parameterization scheme, which tracks the grounding line position within the element (Seroussi et al., 2014).

Model calculations are performed on a horizontally unstructured grid with a higher resolution, l_min=1 km, in fast flow regions and with a coarser resolution, l_max=20 km, in the interior. The vertical discretization comprises 15 layers refined towards the base where vertical shearing becomes more important. The complete mesh comprises 574 056 elements. Velocity, enthalpy (i.e. temperature and microscopic water content), and geometry fields are computed on each vertex of the mesh using piecewise-linear finite elements. The Courant–Friedrichs–Lewy condition (Courant et al., 1928) dictates a time step of 0.025 years. Using the AWI cluster Cray CS400 computer, a simulation with an integration time of 340 years requires ≈8 h on 16 nodes comprised of 36 CPUs.

2.6 Initial state

Future projections of ice sheet evolution first require the determination of the initial state. Different methods are currently used to initialize ice sheets, and it has been shown that the initial state is crucial for projections of ice dynamics (Bindschadler et al., 2013; Goelzer et al., 2018; Nowicki et al., 2013). The recent initMIP–GrIS intercomparison effort (Goelzer et al., 2018) focuses on the different initialization techniques applied in the ice flow modelling community and finds that none of them is the method of choice in terms of a good match to observations and long-term continuity. All methods are required for modelling the projections of the GrIS planned within the CMIP6 phase (Nowicki et al., 2016) on timescales of up to a few hundred years. 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.

Here, we employ a hybrid approach between spin-up and an inversion scheme to estimate the initial state. For the hybrid initialization we make three basic simplifications: (1) the currently observed present-day elevation is taken as constant for the entire glacial cycle. (2) The basal friction coefficient obtained from the inversion is taken as constant for the past glacial cycle, and (3) the temperature changes from the GRIP record are applied to the whole ice sheet without spatial variations.

The ice sheet geometry (bed, ice thickness, and ice sheet mask) is taken from the mass-conserving BedMachine Greenland dataset v2 (Morlighem et al., 2014). The geometric input for thickness and ice sheet mask is masked to exclude glaciers and ice caps surrounding the ice sheet proper. An initial relaxation run over 50 years assuming no sliding and a constant ice temperature of −20 ^∘C is performed to avoid spurious noise that arises from errors and biases in the datasets. A temperature spin-up is then performed using this time-invariant geometry. As the computationally expensive BP approximation is employed, mesh refinements are made at certain points during the whole initialization procedure (see Table 3). The first mesh sequence starts 125 kyr before the present day and runs up to the year 1960, assuming a spatially constant friction coefficient k²=50 s m⁻¹, and is forced with paleoclimatic conditions. The imposed paleoclimatic conditions consist of a multi-year mean from the years 1960 to 1990 of the RACMO2 product (Ettema et al., 2009) and offset by a spatially constant surface temperature anomaly for the last 125 kyr based on the GRIP surface temperature history derived from the Δ¹⁸O record (Dansgaard et al., 1993). The initial ice temperature at 125 kyr before present is a steady-state temperature distribution taken from a spin-up with time-independent climatic conditions from the reference period 1960–1990. The spin-up is done up to 1960 to start the projections before the critical time of Greenland's peripheral glacier mass balance decrease (Noël et al., 2017) with an additional buffer of approximately 30 years.

In the subsequent basal friction inversion, the ice rheology is kept constant using the enthalpy field from the end of the temperature spin-up. The inversion approach infers the basal friction coefficient k² in Eq. (6) by minimizing a cost function that measures the misfit between observed and modelled horizontal velocities (Morlighem et al., 2010). Observed horizontal surface velocities are taken from Rignot and Mouginot (2012). 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 functions are taken from Seroussi et al. (2013).

The procedure for temperature spin-up and inversion is repeated on the subsequent three mesh sequences. The repeated temperature spin-ups start 125, 25, and 15 kyr before 1990 and are again run up to the year 1960. The initial values for the temperature field at these times are taken from the respective times from the previous mesh sequence; the basal friction coefficient is updated from the inversion on the previous mesh sequence. The mesh sequencing reduces the expense of initialization and produces a sufficiently consistent result in terms of velocity and enthalpy. Note that mesh sequences 1–3 are only used during initialization, while the final solution of mesh sequence 4 at year 1960 of this procedure is used as the initial state for all projections presented below.

Table 3 Mesh statistics.

Please note that similar results from this procedure have been submitted to the ISMIP6 initMIP Greenland effort (Goelzer et al., 2018), but the simulations were run with the geothermal flux distribution by Shapiro and Ritzwoller (2004) and additionally with a time-independent climate forcing representing present-day conditions. However, by using the modified heat flux distribution by Greve (2005), we found a generally better agreement with measured basal temperatures at ice core locations. Basically, the comparison of simulated to observed temperatures at the ice base shows too-low temperatures for some locations. As the applied inversion technique for the friction coefficient allows sliding everywhere, the portion of deformational shearing may be underestimated; this cannot be proven without any observations of basal velocities, which unfortunately do not exist. However, for our projections on centennial timescales this is a negligible effect (Seroussi et al., 2013).

2.7 Synthetic and dynamic surface mass balance parameterization

As we perform a one-way coupling of the climatic forcing, the SMB–elevation feedback needs to be considered. Here we rely on the dynamic SMB parameterization developed by Edwards et al. (2014a, b) and previously applied by Goelzer et al. (2013). This relationship was estimated from a set of MAR simulations in which the ice sheet surface elevation was altered. The parameterization assumes that the effect of SMB trends follows a linear relationship:

$$\begin{array}{ll}{\displaystyle }{\mathrm{SMB}}_{\mathrm{dyn}}(x,y,t)& {\displaystyle }={\mathrm{SMB}}_{\mathrm{clim}}(x,y,t)\\ \text{(7)}& {\displaystyle }& {\displaystyle }+b_i\left(h_{\mathrm{s}}(x,y,t)-h_{\mathrm{fix}}(x,y)\right),\end{array}$$

where SMB${}_{\mathrm{dyn}}(x,y,t)$ and SMB${}_{\mathrm{fix}}(x,y,t)$ are the SMB values with and without taking height changes into account, respectively. The surface elevation changes are taken from the ISSM elevation $h_{\mathrm{s}}(x,y,t)$ while running the simulation and a reference elevation h_fix(x, y). In our set-up the reference elevation corresponds to the ISSM ice surface elevation at the initial state.

In this parameterization the SMB gradient b_i is dependent on both location and sign. It can take four values and a separation is made on the location relative to 77^∘ N and on the sign of the SMB. This separates regions of largely different sensitivity, namely the ablation zone with a larger gradient compared to the accumulation zone and a more sensitive ablation zone in the south compared to the north. While a complete uncertainty analysis is given by Edwards et al. (2014a), only the maximum likelihood gradient set, $\mathit{b}=(b_p^N,b_n^N,b_p^S,b_n^S)$, is used here:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{b}_p^N=\mathrm{0.085}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{a}}^{-\mathrm{1}},\\ {\displaystyle }& {\displaystyle }{b}_n^N=\mathrm{0.543}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{a}}^{-\mathrm{1}},\\ {\displaystyle }& {\displaystyle }{b}_p^S=\mathrm{0.063}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{a}}^{-\mathrm{1}},\\ {\displaystyle }& {\displaystyle }{b}_n^S=\mathrm{1.890}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{a}}^{-\mathrm{1}},\end{array}$$

where the subscripts (p, n) and the superscripts (N, S) indicate the evaluation of the SMB sign and the region separation, respectively. Please note that the employed relationships with their parameters may change using a set-up from SEMIC.

A shortcoming of the performed hybrid initialization is that usually a fixed initial ice sheet causes a model drift when imposing the ice thickness equation. This is a result of using an ice sheet that is not in equilibrium with the applied SMB and ice flux divergence. We utilize the local ice thickness imbalance once the ice sheet is released from its fixed topography from a single year unforced relaxation run, i.e. $\mathrm{\Delta }\mathrm{SMB}(x,y,t)=\mathrm{0}$ in Eq. (5). The resulting $\partial H/\partial t$ is subtracted as a surface mass balance correction, SMB_corr(x, y), for all further runs (similar to Goelzer et al., 2018; Price et al., 2011). However, instead of assuming a zero SMB anomaly, one could calculate the anomaly with a GCM input from the CMIP5 pre-industrial scenario. But given the small temperature changes, the SMB anomaly will be close to zero and the calculated ice thickness imbalance is unlikely to be affected by it. However, the final SMB correction is on average 0.01 m a⁻¹, with 5 % of the total ice sheet area having a correction of >25 m a⁻¹, predominantly at marine-terminated ice margins and ice streams (Fig. 6). For these locations, the synthetic SMB correction can be considered additional ice thinning or thickening from dynamic discharge that is not intrinsically simulated. A performed control run with the imposed SMB correction exhibits a small model drift in terms of sea-level equivalent (SLE; black dashed line in Fig. 11 and Sect. 3.3).

The final surface mass balance that the numerical ice flow model sees is composed of several components:

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }\mathrm{SMB}={\mathrm{SMB}}_{\mathrm{clim}}(x,y,t)-{\mathrm{SMB}}_{\mathrm{corr}}(x,y)+{\mathrm{SMB}}_{\mathrm{dyn}}(x,y,t).\end{array}$$

Figure 6 Synthetic surface mass balance SMB_corr calculated from a single year unforced relaxation run (truncated at −25 and 25 m a⁻¹). As the SMB_corr will be subtracted in Eq. (9) positive values represent enforced thinning, and negative values thickening.

3 Results

3.1 Forcing fields

For the different GCMs used we compute ice surface temperature T_s differences between 2100/2300 and 2000 as a multi-year mean over 5 years to reduce the inter-annual variability (Fig. 7). HadGEM2 leads to an increase in temperatures along the northern margins of up to 4 ^∘C. By 2100 the western areas and the vast majority of the ice sheet exceed 2 ^∘C of warming. The only pronounced warming by 2300 is in the north-western regions, while the ice sheet surface temperatures decrease compared to 2100. IPSL exhibits a significantly different pattern with pronounced warming in the centre (up to 3 ^∘C) and in the south-east (up to 4 ^∘C) of the ice sheet. The northern areas reveal moderate warming around 1 ^∘C by 2100. The pattern is similar in 2300, with a moderate cooling in the west compared to 2100. The least warming is found in MIROC5, which even exhibits cooling in the southern areas by about −1 ^∘C in 2100; warming of +1 ^∘C is only reached in the north. By 2300 the entire ice sheet experiences warming; however, this warming is quite moderate compared to the other two GCMs. The low magnitude of warming over Greenland compared to global warming let us infer that the mechanisms of Arctic amplification are not well represented in MIROC5.

Figure 7 Comparison of multi-year mean surface temperature (T_s) differences between 2100–2000 (a–c) and 2300–2000 (d–f) for (a, d) SEMIC-HadGEM2, (b, e) SEMIC-IPSL, and (c, f) SEMIC-MIROC5. The black contour line depicts the present-day ice mask.

Although we do not have a measure to judge future climate warming trends, with respect to the Arctic amplification phenomena the most plausible distribution and magnitude of surface warming is produced by HadGEM2. By contrast, MIROC5 produces less-pronounced warming over Greenland that is similar to the global mean warming but exhibits a plausible pattern of warming. IPSL is spatially and temporally experiencing the largest warming; however, the distribution is not in agreement with the Arctic amplification. Still, the assessment of the GCMs is in line with skill tests performed by Watterson et al. (2014) on a global scale. They assigned skill scores by comparing individual GCM output data against reanalysis data. The analysis indicates that all 25 models have a substantial degree of skill; however, HadGEM2 is ranked in the top, MIROC5 in the middle, and IPSL in the lower part.

Figure 8 presents, in a similar fashion as Fig. 7, the differences in SMB between 2100/2300 and 2000 as a multi-year mean over 5 years each. The difference in SMB 2100–2000 of SEMIC-HadGEM2 indicates a similar pattern as presented by Krapp et al. (2017) using MAR (Fettweis et al., 2013): increasing SMB in the eastern part of the ice sheet with a maximum in the southern half of the ice sheet; at the ice sheet margins ablation is increased. The same pattern is characteristic for 2300–2000, but with a slight decrease in melting and accumulation. The SMB is reduced in the centre, leaving a wide area with differences in SMB of 0.5 m a⁻¹ and less. The SMB difference of SEMIC-IPSL shows a similar pattern with enhanced amplitudes compared to SEMIC-HadGEM2, in particular, and the south-western margin; melting in the south-west is increased by up to 1 m a⁻¹. In contrast, an SMB gain is concentrated in the centre-east by 2300. The most astonishing result is the ΔSMB pattern in SEMIC-MIROC5: increasing the SMB along the south-western and southern margins in contrast to gently decreasing the SMB in the centre of the ice sheet. By 2300 the ΔSMB pattern changes slightly and the SMB is decreasing in the south-western margins. The magnitude of ΔSMB is lower compared to SEMIC-HadGEM2 and SEMIC-IPSL.

Figure 8 Comparison of multi-year mean surface mass balance (SMB) differences between 2100–2000 (a–c) and 2300–2000 (d–f) for (a, d) SEMIC-HadGEM2, (b, e) SEMIC-IPSL, and (c, f) SEMIC-MIROC5. The black contour line depicts the present-day ice mask.

Figure 9 Present-day velocities (year 2000) using SEMIC-HadGEM2: (a) observed velocities, (b) simulated velocities. Observed velocities: Rignot and Mouginot (2012).

Figure 10 Scatter plots of the present-day state (year 2000) using the SMB forcing SEMIC-HadGEM2: (a) velocities, (b) ice surface elevation. Blue and red dots in (a) represent floating and grounded points, respectively. Observed velocities: Rignot and Mouginot (2012); observed surface elevation: Morlighem et al. (2014). The gray line is the identity line.

Figure 11 Sea-level equivalent (SLE, in millimetres) until the year 2100 (left panel) and 2300 (right panel) under RCP2.6 forcing (solid lines) and RCP2.6 forcing without overshoot (dotted-dashed). Additionally, the control run (black dashed line) and the model mean and RMS deviation from Fürst et al. (2015) are shown. The coloured dots represent the onset years of overshooting 1.5 ^∘C in the global mean near-surface air temperature in a 30-year moving window relative to pre-industrial levels.

3.2 Present-day elevation and velocities

Figure 9 displays, as an example, the observed and simulated velocities for the year 2000 (defined here as present day) after a period of forcing with SEMIC-HadGEM2 from 1960 onwards. The resulting horizontal velocity field captures all major features well, including the NEGIS. Outlet glaciers terminating in narrow fjords in the south-eastern region are resolved; however, slow-moving areas tend to retreat below minimum ice thickness and with that the ice extent in this area is underestimated. However, ice surface elevations agree fairly well (Fig. 10b). In general large outlet glaciers like Kangerlussuaq, Helheim, and Jakobshavn Isbræ reveal lower velocities in their fast termini that reflect the high RMS of about 400 m a⁻¹ (Fig. 10a). The RMS analysis here was done on the native grid with high resolution in fast flow regions and the model was already run 40 years forward in time. Compared to these values, the AWI-ISSM results on the regular 5 km grid given in Goelzer et al. (2018) have a lower RMS value of <20 m a⁻¹.

Fürst et al. (2015)

Table 4 Contribution of the Greenland ice sheet to global sea-level change by 2100 and 2300 in millimetres SLE under the RCP2.6 scenario with and without overshoot.

3.3 Projections of mass change

After passing the assumed critical time of declining SMB of the GrIS and the present-day state, the ice sheet experiences a warming and associated mass loss from a decline in surface mass balance. Projections of the evolution of SLE of the ice sheet under the RCP2.6 scenario until 2100 and 2300 are shown in Fig. 11 for each GCM (solid lines) and Table 4. The simulated volume above floatation is converted into the total amount of global sea-level equivalent (SLE) by assuming an ocean area of about 3.618×10⁸ km². Although the control run shows a small model drift in terms of SLE (−1.4 and −0.7 mm for 2100 and 2300, respectively), the RCP2.6 projected SLE is corrected by the control run. By 2100, the model range of Greenland sea-level contributions is between 21.3 and 38.1 mm with an average of 27.9 mm and by 2300 between 36.2 and 85.1 mm with an average of 53.7 mm. Compared to Fürst et al. (2015) our mean values are lower but still in their model range.

The evolution of the mass change, expressed as sea-level equivalent (Fig. 11), shows distinct behaviours: between 1960 and 2000 there is almost no change for SEMIC-HadGEM2 and SEMIC-IPSL, while SEMIC-MIROC5 gains mass; a change in trend with a minor increase between 2000 and 2015 and a steep increase from then on for SEMIC-HadGEM2 and SEMIC-IPSL; and the SLE increase for SEMIC-MIROC5 is more gentle. The steep rise in SLE for SEMIC-HadGEM2 and SEMIC-IPSL is linked to the steep reduction in SMB for both models at the same time. The kink of SLE in SEMIC-HadGEM2 and SEMIC-IPSL around 2050 is caused by a positive SMB anomaly (compare Fig. 4). SEMIC-MIROC5 also shows this peak in SMB, but slightly later at around 2060. These short-term drops in SLE are linked to positive anomalies in SMB. For SEMIC-HadGEM2 the ice sheet contribution until 2300 generally increases continuously, while for SEMIC-IPSL and SEMIC-MIROC5 the increase levels off. This is an intriguing effect as SEMIC-HadGEM2 and IPSL show a similar behaviour in terms of warming over the GrIS (Fig. 1). In fact, the SMB of SEMIC-IPSL recovers from 2050 onwards (Fig. 4), while the SMB of SEMIC-HadGEM2 remains on a low level.

For the RCP2.6 scenario without overshoot the behaviour of SLE for SEMIC-HadGME2 is similar but with lower values. The SLE for SEMIC-MIROC5 is approximately 5 mm lower by 2100 but approaches the same value at 2300 without attaining a pronounced plateau. A striking feature is the much lower SLE estimated from SEMIC-IPSL, which never exceeds a value of 10 mm and gains mass from about 2225 onwards. The average SLE from all three GCMs is 17.4 mm by 2100 and 37.1 mm by 2300, which is approximately one-third less compared to the RCP2.6 scenario.

Figure 12 Lag (j) of projected sea-level rise per year under RCP2.6 forcing (coloured dots) and the modified RCP2.6 forcing without overshoot (coloured circles) as a mean for a time period similar to the observational period (2002–2014). The solid black line indicates the observed value of 0.73 mm a⁻¹ by Rietbroek et al. (2016) and the dashed line the observed value of 0.40 mm a⁻¹ calculated from RACMO2.3 for the period 2002–2014.

Figure 13 Comparison of multi-year mean surface elevation (h_s) differences under RCP2.6 forcing between 2100–2000 (a–c) and 2300–2000 (d–f) for (a, d) SEMIC-HadGEM2, (b, e) SEMIC-IPSL, and (c, f) SEMIC-MIROC5. The black contour line depicts the present-day ice mask. Positive values represent glacier thinning, and negative values thickening. The data are clipped at an ice thickness of 10 m (gray shaded area).

The observed sea-level contribution between 2002 and 2014 is 0.73 mm a⁻¹ (Rietbroek et al., 2016). In the same period, the simulated contribution is only 0.16 mm a⁻¹ for SEMIC-HadGEM2, 0.17 mm a⁻¹ for SEMIC-IPSL, and the lowest for SEMIC-MIROC5 with 0.13 mm a⁻¹. In order to assess a potential temporal lag between the simulated and observed value, mean values of similar periods are calculated (Fig. 12). None of the models reaches the observed value (solid black line in Fig. 12); HadGEM2 reaches a maximum value of 0.59 mm a⁻¹ 13 years later, SEMIC-IPSL a value of 0.48 mm a⁻¹ 12 years later, and SEMIC-MIROC5 a value of 0.36 mm a⁻¹ 40 years later. For the RCP2.6 scenario without overshoot, the values are smaller. Since a future ocean forcing and calving front retreat is not considered here, the response of the ice sheet is likely underestimated. Comparing the sea-level contributions of each SEMIC-GCM to the sea-level contribution of 0.4 mm a⁻¹ calculated from RACMO2.3 for the same period (dashed black line in Fig. 12) reveals a better agreement. SEMIC-HadGEM2 reaches this value 8 years later for the RCP2.6 scenario with overshoot and 9 years later for the RCP2.6 scenario without overshoot; SEMIC-IPSL reaches this value 10 years later for RCP2.6 with overshoot.

3.4 Ice thickness change and dynamic response

Extensive marginal thinning is experienced by forcing the ice sheet with SEMIC-HadGEM2 and SEMIC-IPSL (Fig. 13). In contrast to the mass loss near the margin, the interior shows thickening; IPSL reveals more thickening in the interior. Generally the large-scale pattern of marginal thinning and central thickening correlates with observations (Helm et al., 2014) except that Petermann and Kangerlussuaq glaciers show an opposite trend. With a forcing of MIROC5 the pattern of the elevation change is different with thinning in the southern centre of the ice sheet; the northern centre experienced thickening. Although thinning occurs at the margin it is less extensive compared to the other GCMs.

The response of ice velocities to RCP2.6 forcing is presented in Fig. 14, in which the change in horizontal surface velocities is shown for all scenarios as a difference between 2100–2000 and 2300–2000 (each as a 5-year mean). For all SEMIC-GCM forcings the ice response shows a fairly similar behaviour. The NEGIS, Hagen Bræ, and the Jakobshavn Isbræ, Helheim, and Ryder glaciers experience acceleration; deceleration is present at Petermann and Kangerlussuaq glaciers. However, the magnitude of response is different across all models. It is most prominent at the western margin where SEMIC-HadGEM2 leads to the strongest acceleration, while SEMIC-MIROC5 leads to the lowest.

Figure 14 Comparison of multi-year mean surface velocity (v) differences under RCP2.6 forcing between 2100–2000 (a–c) and 2300–2000 (d–f) for (a, d) SEMIC-HadGEM2, (b, e) SEMIC-IPSL, and (c, f) SEMIC-MIROC5. The black contour line depicts the present-day ice mask. Positive values represent glacier acceleration, and negative values deceleration. The data are clipped at an ice thickness of 10 m (gray shaded area).

4 Discussion

Fürst et al. (2015) performed a comprehensive ensemble study for a suite of 10 GCMs (HadGEM2-ES, IPSL-CM5A-LR and MIROC5 included) and four different RCP scenarios. For the RCP2.6 scenario they estimate a sea-level contribution of 42.3±18.0 mm by 2100 and 88.2±44.8 mm by 2300. Our averaged result of a sea-level contribution under RCP2.6 forcing is slightly lower but still in their ensemble variability. The resultant projection by Fürst et al. (2015) included contributions from lubrication, marine melt, and SMB coupling, while ours accounts for SMB forcing only. The lubrication effect was diagnosed to have a negligible effect on the overall mass budget, but the oceanic influence on the total ice loss explains about half of the mass loss for RCP2.6. Since a future ocean forcing and calving front retreat is not considered here, the response of the ice sheet is likely underestimated. By 2010 the cumulative ice discharge anomaly for SEMIC-HadGEM2 contributes about 15 % to the ice loss. By 2100 and 2300 the contribution is below 3 % and 7 %, respectively, and becomes negligible. For SEMIC-IPSL and SEMIC-MIROC5 the cumulative effect of ice discharge anomaly is less than 10 % of the total mass budget by 2010 and 2100 but increases towards 17 % by 2300. The different behaviour can be explained by the interaction with the SMB and ice dynamics as the relative importance of outlet glacier dynamics decreases with increasing surface melt (Fürst et al., 2015; Goelzer et al., 2013). Increased ice discharge causes dynamic thinning further upstream, lowering of the ice surface, and thereby the intensification of surface melting due to the associated warming of the near surface. Surface melting in turn competes with the discharge increase by removing ice before it reaches the marine margin. The simulated increase in ice discharge for SEMIC-IPSL and SEMIC-MIROC5 is therefore linked to the recovery of SMB over the course of the 22nd century. Still, the SMB remains the dominant factor for mass loss. The speed-up observed from all scenarios merely transports ice from the interior but is melted before it reaches the ice sheet margin. However, the values for sea-level contribution from this study may serve as a lower bound, as processes (ocean forcing and calving) proven to play a major role in GrIS mass loss are not yet represented by the model.

Additionally, the calculation of the surface mass balance is based on different methods. Fürst et al. (2015) rely on the rather simple and empirically derived PDD scheme, while we use a more advanced energy balance approach. So far the sensitivity of melting to warming of this class of models is not well understood. Comparisons of PDD models and energy balance models have suggested that the former are too sensitive to climate change and produce a larger run-off response (Bougamont et al., 2007; Graversen et al., 2011; van de Wal, 1996). On the other hand, Goelzer et al. (2013) attempted to make a robust comparison and find that a PDD model underestimates sea-level rise by 14 %–31 % compared to MAR. An assessment of the SMB and its impact on sea-level contribution calculated by the PDD scheme in Fürst et al. (2015) and the SEMIC from this study cannot be made because of the strong interaction between ice loss, ice dynamics, and external forcings. As the cumulative discharge rates in the mass budget are higher compared to Fürst et al. (2015), this may indicate a lower SMB forcing. However, compared to other models that participate in the initMIP–GrIS exercise (Goelzer et al., 2018), our set-up is neither on the higher nor the lower spectrum of estimated mass loss. Additionally, we have conducted SeaRISE experiments similar to Bindschadler et al. (2013), which showed us that we are within the spread among the models, in particular for the amplified climatic scenarios C1–C3 (not shown here).

The modified RCP2.6 scenario without overshoot projected a sea-level contribution that is on average about 38 % and 31 % less by 2100 and by 2300, respectively. For SEMIC-HadGEM2 and SEMIC-MIROC5 the partition of the mass budget is relatively similar to the RCP2.6 scenario but with a slightly increased cumulative discharge anomaly. For SEMIC-IPSL the behaviour is more irregular, with the ice sheet gaining mass during the last century as a result of an increasing SMB, which is partly compensated for by enhanced ice discharge of up to 40 %. However, the spread of sea-level contribution is much larger compared to the RCP2.6 scenario. In particular, in 2300 the range of sea-level contribution is between 3.4 and 66.9 mm. The very low estimated contribution of 3.4 mm is a result from the SEMIC-IPSL forcing that predicts a relatively high SMB of 364 Gt yr⁻¹ for the characteristic time of overshooting 1.5 ^∘C (column “1.5 ^∘C reached” in Table 2). The SMB is close to present day and therefore SEMIC-IPSL maintains a geometry close to the present day. In contrast, SEMIC-HadGEM2 has declined to 170 Gt yr⁻¹ and SEMIC-MIROC5 to 288 Gt yr⁻¹. The prolongation of these scenarios was done by repeating the forcing from a time window that reveals a stabilized climate. Repeating the last 30-year forcing field window before the characteristic time is not reasonable because the change in warming is strongest during that period and a stabilized climate would not be reached. In fact, we would generate a non-mitigation pathway scenario with constant warming rates that would have a larger melt and therefore make up a larger part of the sea-level contribution (not shown here).

The generally abated sea-level contribution is in agreement with the inferred threshold in global mean temperature before irreversible ice sheet topography changes occur. The simplified assumption behind this threshold is an integrated SMB over the whole ice sheet that becomes negative (Gregory and Huybrechts, 2006). Fettweis et al. (2013) reported a threshold of 3.5 ^∘C relative to pre-industrial levels, which is never exceeded under the RCP2.6 scenario. Assuming a steady-state ice sheet SMB of 400 Gt yr⁻¹ within the reference period, the decline in SMB must be larger than −400 Gt yr⁻¹ to get a continuous, retreating ice sheet margin. If the mean SMB of the GrIS remains positive, a new steady-state ice sheet geometry may be possible, but would require a balancing with the ice outflow.

At last we want to discuss whether studying RCP2.6 allows us to draw significant conclusions on the development of sea-level rise due to mass loss in Greenland. We found that only a fraction of the current observed mass loss in the first 2 decades is represented by the model in RCP2.6. This can be attributed to different factors: first, the current emissions are above the RCP2.6 limit and hence the natural system evolves on a different route than RCP2.6. Secondly, the three GCMs are quite different in response to the RCP2.6 forcing, and the ISM itself does not represent all mechanisms; in particular, the lack of oceanic forcing is causing a reduced sea-level rise. Hence, a new emission scenario that represents the real RCP in the recent past would be most useful for future studies like ours.

5 Conclusions

We have applied climate forcings based on the low emission scenario CMIP5 RCP2.6 of three underlying GCMs (HadGEM2-ES, IPSL-CM5A-LR, MIROC5) to ISSM. Despite all three GCMs being based on RCP2.6, their temperature variation – globally and regionally for the GrIS – is different. Arctic amplification causes a near-surface air temperature increase over Greenland by a factor of ≈2.4 and 2 in HadGEM2-ES and IPSL-CM5A-LR, respectively. MIROC5 reveals nearly no Arctic amplification. In order to force the ice sheet model with a reliable SMB, a physically based surface energy balance model of intermediate complexity (SEMIC) was applied. The estimated sea-level contribution for the RCP2.6 peak and decline scenario in our simulations ranges between 21–38 mm by 2100 and 36–85 mm by 2300 and up to 30 %–40 % higher compared to a scenario without overshoot. Despite the reduced SMB in the warmer climate, a future steady-state ice sheet with lower surface and volume might be possible.

Although the thickness change pattern agrees well with observations and the acceleration of NEGIS, Helheim Glacier, and Jakobshavn Isbræ is captured in our simulations, the estimated sea-level contribution is potentially underestimated due to the following drawbacks of our study: (i) the retreat of glaciers due to oceanic forcing (melt at vertical cliffs and/or calving rates) and (ii) the fact that seasonality due to lubrication arising from supraglacial meltwater is not included. This leads to the conclusion that the projections may serve as a lower bound for the contribution of Greenland to sea-level rise under the RCP2.6 climate scenario. This also limits the advantageous treatment of the physics in our model set-up, meaning that all the benefits from a high-resolution higher-order model are not yet contributing to the extent they potentially could. Our results further indicate that uncertainties stem from the underlying climate model in calculating the surface mass balance.